Analysis of B-cell markers in MCS device patients

Statistical Notebook

Introduction

Raw data

The raw data is shown below. There were 196 rows and 42 columns, consisting of 28 patients sampled over 7 timepoints. However, there is a significant amount of missing data, resulting in only 120 usable datapoints.

missing.ix <- Reduce(intersect, apply(df[,14:42], 2, function(x) which(is.na(x))))
df.raw <- df[-missing.ix,]
df.raw <- df.raw[order(df.raw$PatientID),]
rownames(df.raw) <- 1:nrow(df.raw)
kable(df.raw, 
      digits = 3,
      row.names = T,
      caption = "Raw Data"
) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12) %>%
    scroll_box(width = "100%", height = "300px") 

Raw Data

	PatientID	Age	AgeGreater60	Sex	LowIntermacs	InterMACS	RVAD	Sensitized	VAD Indication	Device Type	Survival	Outcome	Time	num Total PBMC	num lymph	lymph	live lymph	CD3 of live lymph	CD19 of live lymph	CD19+CD27-	CD19+CD27+	CD27+38++plasma blasts	CD27-38++ transitional	CD27-IgD+ mature naive	CD27+IgD- switched memory	CD27-IgD- switched memory	CD27+IgD+ unswitched memory	CD27+IgD-IgM+ switched memory	CD27+IgD+IgM+ nonswitched memory	CD19+27+IgG+IgM- memory	CD19+24dim38dim naive mature	CD19+24+38++transitional	CD19CD24hiCD38-memory	CD19+27-38+CD5+transitionals	CD19+CD268+	CD268 of +27-38++transitional	CD19+CD11b+	CD19+CD5+	CD19+CD27+CD24hi	CD19+CD5+CD24hi	CD19+CD5+CD11b+	CD19+27+IgD-38++IgG ASC
1	1	65	older	Male	High	2	No	NA	BTT	HMII	alive	Alive s/p OHT	0	169154	35496	20.98	99.53	25.37	19.82	86.02	13.98	1.26	2.54	57.81	11.74	28.06	2.39	21.74	14.36	0.72	81.05	2.86	13.60	0.58	95.22	84.83	10.47	4.33	8.50	1.39	3.13	1.09
2	1	65	older	Male	High	2	No	NA	BTT	HMII	alive	Alive s/p OHT	1	63082	9915	15.72	99.26	32.45	26.26	90.06	9.94	0.89	2.09	54.27	7.70	35.71	2.32	31.10	20.47	1.57	84.06	2.59	10.79	0.26	97.02	92.59	7.50	3.13	6.42	1.32	2.44	1.49
3	1	65	older	Male	High	2	No	NA	BTT	HMII	alive	Alive s/p OHT	3	75921	21721	28.61	99.31	24.86	33.01	86.38	13.62	1.54	1.26	45.74	10.98	40.44	2.84	31.80	18.04	1.15	86.21	4.82	6.56	0.24	90.41	73.33	10.00	7.25	10.12	3.22	5.67	2.05
4	1	65	older	Male	High	2	No	NA	BTT	HMII	alive	Alive s/p OHT	8	213808	36002	16.84	99.19	35.95	14.55	67.94	32.06	4.06	2.17	27.92	28.92	39.64	3.52	20.31	9.10	1.33	68.69	7.10	20.42	1.88	87.38	44.25	10.78	11.12	19.15	5.54	7.35	1.07
5	1	65	older	Male	High	2	No	NA	BTT	HMII	alive	Alive s/p OHT	21	106970	31867	29.79	99.18	42.36	15.02	84.69	15.31	1.94	2.11	55.41	13.14	29.13	2.32	20.22	12.88	2.08	83.21	5.20	9.50	1.10	94.44	73.00	9.04	6.99	9.22	2.80	4.78	1.44
6	2	81	older	Male	Low	3	No	NA	DT	HMII	dead	Died	0	119377	77509	64.93	99.61	58.50	7.99	62.67	37.33	9.10	4.82	36.73	31.90	25.62	5.74	11.47	8.95	0.40	69.54	7.69	9.41	2.61	76.44	57.91	23.29	20.76	15.39	5.42	14.90	2.60
7	2	81	older	Male	Low	3	No	NA	DT	HMII	dead	Died	3	96940	57658	59.48	99.70	55.52	9.21	72.15	27.85	3.85	3.82	48.41	21.75	23.49	6.35	12.10	16.04	0.28	75.30	3.72	6.65	2.11	80.52	67.33	15.80	15.23	6.71	1.32	8.86	0.26
8	2	81	older	Male	Low	3	No	NA	DT	HMII	dead	Died	5	154373	73969	47.92	99.65	45.18	8.77	67.36	32.64	6.59	2.77	39.10	26.28	27.95	6.67	17.89	15.71	0.39	71.57	6.26	9.88	1.78	81.86	72.07	20.12	20.17	14.80	5.99	14.65	3.71
9	2	81	older	Male	Low	3	No	NA	DT	HMII	dead	Died	8	155610	63375	40.73	99.51	38.42	12.40	71.69	28.31	4.36	2.62	37.82	23.43	33.50	5.24	17.37	14.92	0.60	65.04	2.80	11.23	2.54	83.55	37.56	15.55	15.45	9.55	2.66	10.91	0.28
10	2	81	older	Male	Low	3	No	NA	DT	HMII	dead	Died	14	244245	115109	47.13	99.47	47.40	9.37	67.50	32.50	3.99	2.08	50.08	26.01	17.22	6.68	15.34	15.25	0.23	77.86	3.21	10.49	1.48	82.13	59.19	21.45	19.77	14.45	2.57	14.94	0.94
11	3	58	younger	Male	Low	3	Yes	No	BTT	HMII	alive	Alive s/p OHT	1	246811	86529	35.06	93.49	67.96	17.48	81.66	18.34	2.05	6.85	51.33	13.69	29.91	5.08	15.68	5.94	0.19	61.95	0.76	13.51	4.07	92.62	34.06	6.90	12.07	8.00	2.54	3.85	0.26
12	3	58	younger	Male	Low	3	Yes	No	BTT	HMII	alive	Alive s/p OHT	3	406771	90343	22.21	97.34	49.28	17.02	75.04	24.96	4.47	18.36	28.45	21.86	46.07	3.62	24.97	4.33	0.13	47.41	1.64	19.70	18.15	72.76	8.08	25.69	28.35	13.95	4.66	19.20	0.41
13	3	58	younger	Male	Low	3	Yes	No	BTT	HMII	alive	Alive s/p OHT	21	330191	81636	24.72	94.57	45.14	14.70	70.51	29.49	8.94	16.23	38.92	23.53	30.83	6.72	23.93	7.07	0.18	46.63	2.35	16.12	8.59	74.45	39.90	24.19	21.18	12.62	3.26	12.93	0.38
14	4	67	older	Male	Low	3	No	Yes	BTT	HVAD	alive	Alive s/p OHT	0	309890	41896	13.52	97.03	47.13	8.08	77.31	22.69	2.74	3.29	59.03	19.80	19.37	1.80	21.10	4.46	19.79	71.98	2.41	17.39	0.20	95.80	75.00	5.79	1.89	16.45	0.91	0.67	0.77
15	4	67	older	Male	Low	3	No	Yes	BTT	HVAD	alive	Alive s/p OHT	1	274228	15653	5.71	97.85	35.76	28.87	51.45	48.55	2.60	2.13	30.85	43.13	22.77	3.26	20.09	3.92	30.51	57.58	1.31	28.20	0.18	91.23	59.57	12.08	2.04	34.44	0.38	0.16	0.16
16	4	67	older	Male	Low	3	No	Yes	BTT	HVAD	alive	Alive s/p OHT	3	591939	57589	9.73	96.76	38.20	42.54	75.71	24.29	2.14	1.33	53.03	21.26	24.16	1.55	25.50	4.02	20.27	67.18	1.14	21.71	0.11	95.28	66.98	3.79	2.43	17.70	0.65	0.84	0.18
17	4	67	older	Male	Low	3	No	Yes	BTT	HVAD	alive	Alive s/p OHT	5	189634	46711	24.63	97.74	34.84	37.32	75.30	24.70	0.70	0.55	46.97	22.05	30.28	0.70	23.67	1.56	31.44	73.33	0.75	18.59	0.00	84.61	86.17	3.44	2.11	18.68	0.94	0.63	0.69
18	4	67	older	Male	Low	3	No	Yes	BTT	HVAD	alive	Alive s/p OHT	14	182319	40475	22.20	96.28	43.36	11.77	69.34	30.66	1.59	0.89	52.64	27.32	18.56	1.48	22.79	2.88	28.61	75.47	0.33	14.54	0.06	95.90	53.66	7.81	3.27	21.04	1.57	1.50	0.72
19	5	65	older	Female	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	0	1073919	336538	31.34	100.00	63.29	14.80	93.05	6.95	0.73	8.05	75.49	4.55	19.16	0.80	7.25	8.78	28.04	88.72	2.51	5.18	1.02	90.49	68.18	2.32	4.09	1.69	3.48	1.27	5.78
20	5	65	older	Female	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	1	270301	97229	35.97	91.63	18.43	61.67	97.68	2.32	0.20	5.85	69.45	1.63	28.19	0.73	27.57	21.43	17.05	84.32	7.26	2.73	0.84	99.44	97.60	0.54	2.64	1.51	3.63	0.29	3.40
21	5	65	older	Female	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	3	499016	78019	15.63	89.41	25.83	51.38	92.98	7.02	0.32	4.66	51.34	6.31	41.59	0.76	36.49	5.96	22.81	73.48	6.73	7.93	0.33	98.78	96.52	1.31	2.06	5.19	2.68	0.59	2.45
22	5	65	older	Female	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	8	115825	32513	28.07	91.04	44.61	15.76	80.52	19.48	4.16	7.82	49.91	18.37	30.52	1.20	14.87	3.33	30.19	61.08	8.70	12.56	1.58	89.80	76.71	2.87	6.07	11.12	8.08	1.59	1.17
23	5	65	older	Female	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	21	407080	97651	23.99	91.29	66.96	12.20	84.66	15.34	2.88	5.60	61.89	13.73	22.65	1.74	20.16	7.23	29.39	68.49	5.28	10.60	1.17	92.50	54.35	3.48	6.53	9.61	8.34	1.85	2.83
24	6	43	younger	Male	High	2	Yes	No	BTT	HMII	alive	Alive s/p OHT	0	515760	266115	51.60	96.30	69.66	8.57	72.22	27.78	6.26	4.77	43.09	26.01	28.26	2.63	16.67	2.67	28.25	50.35	0.41	12.78	4.47	84.76	19.47	9.56	4.55	8.29	0.87	2.18	0.75
25	6	43	younger	Male	High	2	Yes	No	BTT	HMII	alive	Alive s/p OHT	8	539521	223872	41.49	96.72	60.02	12.63	80.13	19.87	1.66	2.48	43.56	19.52	35.24	1.68	20.46	3.67	20.54	55.76	0.18	10.80	2.43	87.82	13.25	6.55	5.96	4.65	0.95	2.51	0.65
26	6	43	younger	Male	High	2	Yes	No	BTT	HMII	alive	Alive s/p OHT	21	564317	216312	38.33	97.40	55.80	4.83	74.19	25.81	5.35	5.77	37.90	25.72	34.76	1.62	17.99	2.18	31.75	61.61	0.77	13.44	6.97	76.55	6.97	10.34	11.14	4.46	0.94	5.14	2.43
27	7	49	younger	Male	Low	3	No	No	BTT	HVAD	alive	Alive s/p OHT	0	89251	65181	73.03	76.76	84.53	1.15	74.65	25.35	1.91	1.39	55.03	22.40	19.62	2.95	15.65	6.80	25.85	70.49	0.87	23.26	0.23	83.16	25.00	8.85	5.03	11.63	1.56	2.95	3.08
28	7	49	younger	Male	Low	3	No	No	BTT	HVAD	alive	Alive s/p OHT	3	375167	234817	62.59	98.46	60.29	10.18	87.72	12.28	0.46	0.48	68.54	9.90	19.03	2.54	30.58	19.69	15.54	72.77	0.13	20.35	0.03	92.56	36.28	3.75	2.70	9.17	1.19	1.27	0.63
29	7	49	younger	Male	Low	3	No	No	BTT	HVAD	alive	Alive s/p OHT	5	256306	162176	63.27	93.82	68.43	8.45	85.45	14.55	0.70	0.53	71.99	11.98	13.37	2.66	32.41	17.96	12.46	69.74	0.41	24.63	0.02	95.72	44.12	4.47	2.23	12.99	1.30	1.20	0.77
30	7	49	younger	Male	Low	3	No	No	BTT	HVAD	alive	Alive s/p OHT	8	463136	330714	71.41	98.79	70.51	7.30	87.51	12.49	1.78	0.58	69.60	10.61	17.82	1.97	29.83	14.55	12.59	75.03	0.15	21.61	0.07	93.59	30.94	4.05	3.07	10.00	1.08	1.35	0.62
31	8	43	younger	Male	Low	3	No	No	BTT	HMII	alive	Alive s/p OHT	0	92221	27988	30.35	85.45	55.78	1.51	85.56	14.44	0.28	0.28	55.28	4.72	39.72	0.28	0.00	1.89	1.89	89.44	0.00	2.78	0.00	15.00	0.00	3.33	0.83	0.28	0.28	3.33	0.00
32	8	43	younger	Male	Low	3	No	No	BTT	HMII	alive	Alive s/p OHT	1	273746	46755	17.08	93.25	50.14	7.57	88.27	11.73	0.09	2.12	54.52	3.36	42.00	0.12	8.40	6.87	0.00	92.00	0.12	3.45	0.00	3.94	11.43	1.12	0.64	2.76	0.33	2.15	0.00
33	8	43	younger	Male	Low	3	No	No	BTT	HMII	alive	Alive s/p OHT	3	156053	71752	45.98	78.38	35.43	27.12	78.57	21.43	0.24	0.60	58.56	5.80	18.33	17.31	12.67	72.96	0.53	94.06	0.49	5.11	0.28	22.79	27.47	0.81	5.53	0.01	0.29	0.19	0.00
34	8	43	younger	Male	Low	3	No	No	BTT	HMII	alive	Alive s/p OHT	8	156053	71752	45.98	78.38	35.43	27.12	78.57	21.43	0.24	0.60	58.56	5.80	18.33	17.31	12.67	72.96	0.53	94.06	0.49	5.11	0.28	22.79	27.47	0.81	5.53	0.01	0.29	0.19	0.00
35	9	60	younger	Female	High	2	Yes	Yes	BTT	PVAD	dead	Died	3	116462	41076	35.27	95.44	33.79	25.24	84.51	15.49	2.92	0.99	45.34	13.19	39.12	2.36	17.56	8.75	19.73	79.35	5.17	10.80	0.71	90.04	29.59	10.16	6.81	8.21	3.37	4.86	4.47
36	9	60	younger	Female	High	2	Yes	Yes	BTT	PVAD	dead	Died	5	222650	85649	38.47	98.40	32.09	19.82	77.54	22.46	5.02	1.10	43.96	19.20	33.50	3.34	17.47	9.77	15.63	71.89	6.68	14.08	0.70	86.49	27.87	13.04	9.84	12.30	4.35	7.98	6.39
37	9	60	younger	Female	High	2	Yes	Yes	BTT	PVAD	dead	Died	8	102274	63057	61.65	97.73	27.89	20.66	73.75	26.25	7.33	3.45	35.79	21.92	37.73	4.56	9.97	8.20	16.96	63.96	2.40	11.16	3.07	73.97	5.92	16.69	17.77	11.74	4.90	12.19	6.09
38	9	60	younger	Female	High	2	Yes	Yes	BTT	PVAD	dead	Died	14	133530	81665	61.16	96.46	22.03	16.15	35.30	64.70	13.74	2.33	8.48	52.28	26.37	12.87	5.29	5.79	19.15	47.48	5.65	17.83	5.43	39.88	6.08	52.08	62.11	37.45	21.95	50.15	3.84
39	10	36	younger	Male	High	2	Yes	NA	BTT	PVAD	alive	Alive s/p OHT	0	481928	162052	33.63	99.70	64.73	5.05	3.91	96.09	8.89	1.34	0.72	94.99	3.14	1.15	4.07	0.84	0.32	90.28	2.84	1.58	1.27	6.82	9.17	14.61	90.11	4.95	5.71	13.01	0.51
40	10	36	younger	Male	High	2	Yes	NA	BTT	PVAD	alive	Alive s/p OHT	8	38113	21486	56.37	97.86	58.44	7.88	20.47	79.53	8.09	1.69	1.21	75.00	18.96	4.83	10.07	2.59	0.53	82.19	2.29	9.00	3.59	15.16	0.00	23.07	63.95	4.77	6.52	16.91	0.40
41	10	36	younger	Male	High	2	Yes	NA	BTT	PVAD	alive	Alive s/p OHT	21	80991	19210	23.72	99.00	50.61	3.50	12.33	87.67	14.74	4.81	1.05	85.41	11.13	2.41	8.48	0.52	1.56	80.00	8.27	3.76	7.41	13.68	0.00	30.83	77.74	10.83	14.14	19.85	2.81
42	11	61	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	0	411616	119210	28.96	90.82	73.63	5.88	74.30	25.70	2.80	6.46	57.94	19.20	18.72	4.13	19.10	15.42	15.30	53.51	13.78	25.19	0.21	47.23	53.77	5.75	5.92	22.63	3.46	2.55	0.74
43	11	61	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	1	430509	79492	18.46	89.31	52.10	20.45	87.13	12.87	1.23	4.76	58.11	9.15	30.55	2.20	29.86	18.05	10.37	63.63	10.23	17.43	0.06	74.05	72.50	2.36	2.98	10.68	1.72	1.03	2.89
44	11	61	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	3	395713	119670	30.24	97.90	23.79	6.27	82.17	17.83	0.79	1.99	49.95	13.01	33.91	3.13	28.36	16.83	5.00	56.33	7.13	26.73	0.02	79.50	82.88	1.96	5.55	16.30	3.61	1.88	0.42
45	11	61	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	8	329840	73468	22.27	98.63	70.61	3.73	51.50	48.50	8.32	3.77	33.46	41.29	20.92	4.33	25.19	5.89	24.07	46.40	3.73	34.64	0.07	64.14	28.43	10.35	9.61	38.78	4.62	5.55	2.91
46	11	61	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	14	350294	88823	25.36	97.63	74.37	4.37	73.35	26.65	4.91	4.83	56.54	21.32	19.84	2.30	34.68	11.25	18.37	61.98	5.78	20.55	0.04	63.91	47.54	6.81	6.17	21.77	2.90	2.74	3.00
47	12	38	younger	Male	High	2	Yes	Yes	BTT	TAH	dead	Died s/p OHT	0	128451	60927	47.43	96.82	60.77	4.70	77.26	22.74	2.17	1.95	55.58	18.40	21.54	4.47	27.86	18.89	15.33	69.90	0.94	24.40	0.14	70.73	1.85	8.01	8.66	20.21	5.45	3.32	0.39
48	12	38	younger	Male	High	2	Yes	Yes	BTT	TAH	dead	Died s/p OHT	1	142796	22573	15.81	80.67	47.31	2.28	63.37	36.63	9.16	6.51	26.02	33.49	37.11	3.37	38.13	7.50	13.13	53.25	5.78	26.27	1.57	40.24	0.00	27.71	22.17	29.16	18.80	15.66	8.11
49	12	38	younger	Male	High	2	Yes	Yes	BTT	TAH	dead	Died s/p OHT	3	157167	38008	24.18	93.47	34.86	3.95	69.71	30.29	3.92	3.92	33.64	26.30	35.78	4.28	39.54	11.49	14.94	66.86	2.85	24.31	1.04	64.50	7.27	21.60	20.53	27.66	14.54	14.18	2.12
50	12	38	younger	Male	High	2	Yes	Yes	BTT	TAH	dead	Died s/p OHT	5	121823	23097	18.96	93.48	2.51	2.87	99.52	0.48	0.00	0.00	0.97	0.48	98.55	0.00	50.00	25.00	25.00	94.35	0.00	2.58	0.00	0.16	0.00	17.45	87.40	0.00	0.16	18.26	0.00
51	12	38	younger	Male	High	2	Yes	Yes	BTT	TAH	dead	Died s/p OHT	8	196515	88579	45.07	85.54	17.17	3.28	59.38	40.62	3.21	1.41	35.44	35.28	23.87	5.42	62.85	12.32	3.71	66.05	4.38	23.22	0.27	74.53	31.43	25.51	28.12	39.78	25.79	21.57	0.79
52	13	61	older	Male	High	1	Yes	No	BTT	TAH	alive	Alive s/p OHT	1	87827	15751	17.93	76.38	31.60	12.17	77.46	22.54	6.35	2.53	50.34	20.22	27.05	2.39	19.68	7.62	19.05	75.55	7.38	9.49	0.71	94.26	89.19	11.68	17.96	18.78	6.97	10.86	12.46
53	13	61	older	Male	High	1	Yes	No	BTT	TAH	alive	Alive s/p OHT	3	55539	13608	24.50	79.00	28.02	15.78	81.90	18.10	6.13	2.42	66.92	14.33	14.98	3.77	20.81	17.45	18.12	78.01	7.08	8.61	0.72	94.04	70.73	11.32	16.69	15.21	5.78	10.50	16.39
54	13	61	older	Male	High	1	Yes	No	BTT	TAH	alive	Alive s/p OHT	8	295005	36785	12.47	94.58	36.55	9.14	86.32	13.68	2.45	1.82	37.01	12.52	49.18	1.29	16.24	4.71	12.47	8.36	1.73	6.64	0.81	84.78	46.55	12.52	9.28	5.94	1.73	4.47	2.49
55	13	61	older	Male	High	1	Yes	No	BTT	TAH	alive	Alive s/p OHT	14	60695	20471	33.73	94.91	6.93	7.42	83.97	16.03	3.40	0.69	24.36	14.71	59.33	1.60	18.43	5.07	11.52	54.82	2.91	16.10	0.50	82.44	30.00	51.08	13.25	10.41	3.68	11.31	6.06
56	14	39	younger	Female	Low	3	Yes	No	BTT	HMII	dead	Died post OHT	0	52016	30920	59.44	95.77	55.57	24.62	58.73	41.27	2.14	0.97	53.01	5.05	5.45	36.50	9.39	81.29	2.88	81.48	5.65	9.02	0.48	97.75	71.83	0.67	22.89	0.34	1.78	0.49	NA
57	14	39	younger	Female	Low	3	Yes	No	BTT	HMII	dead	Died post OHT	1	75516	48801	64.62	92.13	52.90	18.64	70.39	29.61	2.92	2.21	62.18	5.13	7.98	24.70	13.79	69.54	3.19	84.14	3.64	6.36	0.41	94.93	43.78	0.57	15.61	0.23	0.91	0.36	NA
58	14	39	younger	Female	Low	3	Yes	No	BTT	HMII	dead	Died post OHT	3	250612	130454	52.05	95.75	60.05	3.29	46.04	53.96	1.34	3.89	33.08	32.47	8.29	26.15	3.56	34.61	0.09	76.71	0.39	13.20	0.28	39.79	14.38	2.26	1.56	10.45	2.89	3.04	0.08
59	14	39	younger	Female	Low	3	Yes	No	BTT	HMII	dead	Died post OHT	8	212478	120477	56.70	85.12	48.36	12.51	79.57	20.43	1.35	2.14	58.73	12.87	18.44	9.96	5.63	34.32	0.11	86.74	0.41	6.88	0.15	35.30	53.09	1.23	0.75	2.36	1.04	1.47	NA
60	15	70	older	Female	High	2	No	NA	DT	HMII	alive	Alive	0	86135	70306	81.62	93.71	74.50	1.59	61.74	38.26	1.91	2.67	33.97	14.50	26.81	24.71	20.67	62.47	9.03	73.00	2.00	24.71	1.27	20.23	71.43	0.48	12.79	0.38	3.24	0.48	4.60
61	15	70	older	Female	High	2	No	NA	DT	HMII	alive	Alive	1	81004	41334	51.03	94.06	53.64	2.99	61.70	38.30	3.10	2.75	39.07	13.34	22.12	25.47	12.58	63.11	15.99	78.92	2.93	17.38	3.35	23.24	71.88	0.52	11.70	0.09	1.29	0.52	9.83
62	15	70	older	Female	High	2	No	NA	DT	HMII	alive	Alive	5	90024	39434	43.80	95.32	68.83	5.62	71.59	28.41	0.85	1.66	35.23	13.02	35.80	15.96	22.94	52.37	5.85	76.42	1.14	21.83	0.62	14.54	48.57	0.38	4.69	0.00	0.80	0.38	0.62
63	15	70	older	Female	High	2	No	NA	DT	HMII	alive	Alive	14	106791	80532	75.41	88.82	79.92	1.57	62.32	37.68	1.70	4.02	28.04	22.86	34.11	15.00	23.92	39.41	15.03	61.61	3.66	34.55	1.97	20.36	73.33	0.80	8.57	0.27	1.79	0.98	0.71
64	16	63	older	Male	High	1	No	NA	BTT	HMII	dead	Died	3	366713	77282	21.07	98.61	66.55	12.44	93.65	6.35	0.53	0.60	62.22	4.68	32.38	0.72	19.94	9.89	28.23	83.18	0.51	14.38	0.30	94.72	54.39	2.07	11.33	2.86	3.40	1.48	3.58
65	16	63	older	Male	High	1	No	NA	BTT	HMII	dead	Died	5	262251	55344	21.10	99.55	62.37	14.39	93.00	7.00	1.03	0.69	65.71	4.74	28.37	1.17	25.93	17.28	25.93	79.87	1.20	17.05	0.48	97.76	58.18	2.27	11.13	4.12	4.91	1.87	6.10
66	16	63	older	Male	High	1	No	NA	BTT	HMII	dead	Died	8	309193	61091	19.76	98.15	56.62	7.24	83.61	16.39	3.27	1.52	58.03	11.90	28.13	1.93	24.65	11.36	1.52	79.53	1.45	13.33	0.97	92.63	25.76	3.50	9.23	7.57	3.73	3.48	0.00
67	17	68	older	Male	Low	3	No	NA	DT	HMII	dead	Died	0	147997	54649	36.93	91.75	79.63	7.62	74.86	25.14	0.60	3.45	62.45	22.87	12.30	2.38	45.98	9.28	13.61	74.33	2.67	21.19	0.84	91.31	90.15	3.58	2.25	21.79	1.36	0.58	0.34
68	17	68	older	Male	Low	3	No	NA	DT	HMII	dead	Died	1	177687	56400	31.74	97.13	69.76	12.60	81.83	18.17	0.65	1.61	55.77	16.13	26.02	2.09	53.19	10.70	11.57	82.92	1.38	13.98	0.28	89.83	85.59	3.78	1.85	14.65	1.30	0.68	0.18
69	17	68	older	Male	Low	3	No	NA	DT	HMII	dead	Died	3	233155	67140	28.80	98.44	74.89	5.32	81.96	18.04	1.00	1.88	49.03	16.59	32.90	1.48	49.92	7.20	12.83	83.44	1.39	12.98	0.14	89.70	69.70	4.58	1.62	13.89	1.08	1.08	0.17
70	17	68	older	Male	Low	3	No	NA	DT	HMII	dead	Died	5	99257	29455	29.68	97.33	62.48	4.73	80.60	19.40	1.84	1.25	45.06	18.36	35.25	1.33	44.61	6.32	12.27	82.89	0.88	13.20	0.09	83.55	35.29	7.23	3.54	14.75	2.29	3.02	1.20
71	17	68	older	Male	Low	3	No	NA	DT	HMII	dead	Died	8	76358	18702	24.49	97.32	29.48	11.10	86.29	13.71	2.43	0.50	37.92	12.13	48.32	1.63	46.45	11.70	11.35	83.96	0.20	10.50	0.06	57.82	50.00	9.95	29.60	10.20	0.64	9.36	0.80
72	18	37	younger	Female	High	2	Yes	NA	BTT	PVAD	dead	Died s/p OHT	0	56374	12141	21.54	67.38	27.12	3.32	51.84	48.16	1.84	1.84	2.94	45.59	47.79	3.68	0.00	0.00	2.29	92.65	0.37	5.51	0.72	30.15	60.00	6.99	9.56	0.37	0.00	3.68	0.00
73	18	37	younger	Female	High	2	Yes	NA	BTT	PVAD	dead	Died s/p OHT	3	7018	3464	49.36	89.03	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
74	18	37	younger	Female	High	2	Yes	NA	BTT	PVAD	dead	Died s/p OHT	8	80210	21828	27.21	85.74	20.38	3.50	61.98	38.02	8.55	2.90	5.50	34.96	56.18	3.36	0.00	0.00	3.21	89.31	0.61	2.60	1.49	23.66	0.00	21.37	9.16	1.07	0.15	5.34	1.31
75	18	37	younger	Female	High	2	Yes	NA	BTT	PVAD	dead	Died s/p OHT	14	55068	19145	34.77	75.16	16.90	3.30	64.00	36.00	4.21	2.11	10.53	34.11	52.42	2.95	0.00	0.00	2.94	92.42	0.00	3.37	1.00	22.74	0.00	12.42	11.16	0.63	0.00	4.63	0.00
76	18	37	younger	Female	High	2	Yes	NA	BTT	PVAD	dead	Died s/p OHT	21	76673	23615	30.80	83.39	26.27	2.06	63.55	36.45	5.42	5.42	5.17	34.48	58.13	2.22	0.00	0.00	2.72	89.41	0.00	3.45	1.95	18.23	4.55	18.23	8.87	1.23	0.25	4.68	0.00
77	19	46	younger	Male	High	2	No	NA	DT	HMII	alive	Alive	0	226762	25868	11.41	89.41	25.31	12.05	55.06	44.94	11.31	2.69	14.79	42.57	39.95	2.69	19.87	3.88	32.88	68.95	10.48	13.57	0.59	80.62	22.67	31.91	19.38	15.58	10.37	16.55	8.18
78	19	46	younger	Male	High	2	No	NA	DT	HMII	alive	Alive	3	244049	39386	16.14	94.94	40.46	13.23	40.93	59.07	13.36	3.17	8.33	54.90	32.28	4.49	18.52	4.21	23.83	55.18	13.87	19.39	1.44	70.53	35.67	33.62	33.47	22.86	15.26	25.81	5.11
79	19	46	younger	Male	High	2	No	NA	DT	HMII	alive	Alive	5	344736	250121	72.55	89.83	7.54	3.81	45.84	54.16	3.96	1.53	13.39	52.70	32.18	1.72	67.14	2.92	2.44	74.68	6.81	13.91	0.21	63.33	58.02	33.05	30.94	19.83	13.96	26.74	0.22
80	19	46	younger	Male	High	2	No	NA	DT	HMII	alive	Alive	8	216934	140093	64.58	96.37	4.78	7.22	22.49	77.51	13.49	0.64	7.22	67.61	15.17	10.00	35.88	7.82	17.22	51.28	39.40	6.56	1.33	62.59	58.06	66.67	64.59	53.00	50.26	63.78	4.97
81	19	46	younger	Male	High	2	No	NA	DT	HMII	alive	Alive	14	215820	66785	30.94	92.54	17.20	6.01	25.88	74.12	9.15	1.45	3.09	71.72	22.44	2.74	65.43	3.42	0.04	69.38	15.93	8.85	0.74	50.58	24.07	52.78	48.88	32.28	28.19	44.04	0.04
82	19	46	younger	Male	High	2	No	NA	DT	HMII	alive	Alive	21	670692	523014	77.98	96.05	5.07	2.64	23.61	76.39	4.11	0.49	4.79	73.53	18.68	3.00	42.40	2.72	4.49	71.95	20.81	5.51	0.45	48.22	44.62	60.13	60.26	39.07	36.82	56.38	1.49
83	20	75	older	Male	High	2	No	NA	DT	HMII	alive	Alive	3	121794	28129	23.10	97.05	52.79	16.42	70.08	29.92	0.29	0.54	50.71	21.46	22.96	4.86	34.18	12.51	23.49	39.78	0.31	55.18	0.35	96.01	16.67	1.74	5.27	23.38	2.95	1.00	0.10
84	20	75	older	Male	High	2	No	NA	DT	HMII	alive	Alive	5	82157	17295	21.05	99.20	72.54	7.41	68.21	31.79	1.10	0.55	48.47	22.90	23.60	5.04	31.90	13.57	18.10	37.84	0.47	56.73	0.35	92.37	14.29	2.05	7.79	26.20	4.72	1.42	1.02
85	20	75	older	Male	High	2	No	NA	DT	HMII	alive	Alive	8	154810	27067	17.48	98.04	72.18	5.03	74.91	25.09	0.90	1.65	58.05	14.76	20.82	6.37	23.43	25.71	22.57	44.94	0.60	49.36	0.71	86.59	22.73	3.30	7.64	20.15	4.87	1.72	0.00
86	20	75	older	Male	High	2	No	NA	DT	HMII	alive	Alive	21	197936	15446	7.80	95.86	50.35	5.44	79.13	20.87	3.98	3.11	52.30	15.90	29.81	1.99	20.23	5.78	24.86	40.62	1.12	50.06	1.89	62.36	8.00	7.20	9.94	13.17	3.60	5.09	1.55
87	21	30	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	0	114586	35996	31.41	100.00	59.03	11.31	75.87	24.13	1.52	7.69	69.83	11.67	10.20	8.30	27.80	29.94	2.04	84.64	0.64	10.54	0.61	91.57	61.34	4.99	7.15	11.57	5.92	2.53	1.92
88	21	30	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	3	349054	129419	37.08	95.25	62.35	8.81	79.64	20.36	2.72	4.45	57.70	13.93	25.89	2.48	57.14	13.49	8.69	74.29	3.76	12.35	1.41	85.87	62.53	3.41	5.43	7.67	1.40	1.92	1.06
89	21	30	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	5	447309	176789	39.52	93.00	53.57	9.32	80.92	19.08	2.16	5.73	56.72	12.51	28.43	2.34	47.68	12.19	11.00	73.80	4.44	12.53	2.57	83.88	49.09	5.24	6.82	8.32	2.18	2.28	1.72
90	21	30	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	8	558290	301468	54.00	95.62	47.73	7.04	79.13	20.87	3.43	7.49	53.19	13.48	30.32	3.01	36.36	14.54	12.10	70.55	5.15	11.89	3.51	79.79	52.89	6.98	8.50	7.82	1.92	3.09	1.96
91	21	30	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	14	479402	362845	75.69	97.95	44.30	6.30	76.53	23.47	3.90	11.10	34.99	12.04	50.31	2.66	21.93	14.88	10.05	70.21	4.66	7.31	7.73	49.96	29.38	13.44	17.37	5.64	3.01	6.95	2.28
92	21	30	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	21	334942	191418	57.15	93.28	70.35	4.04	70.89	29.11	6.20	17.70	56.45	19.57	20.44	3.53	38.91	14.64	10.35	55.11	6.61	14.34	8.70	80.96	58.92	7.58	7.96	10.11	2.08	2.08	2.48
93	22	67	older	Male	High	2	Yes	NA	BTT	HMII	dead	Died	0	129980	41600	32.00	91.40	39.03	2.85	49.63	50.37	7.02	11.55	28.37	29.11	20.79	21.72	19.43	28.34	9.80	69.59	4.90	20.33	6.13	31.15	33.60	20.15	22.27	0.28	4.25	16.54	2.09
94	22	67	older	Male	High	2	Yes	NA	BTT	HMII	dead	Died	1	59478	20182	33.93	89.07	44.50	10.08	60.93	39.07	2.10	2.70	42.88	13.69	17.66	25.77	13.06	46.12	3.13	81.35	1.16	16.11	0.83	37.58	6.12	5.52	11.64	0.06	1.49	2.76	4.07
95	22	67	older	Male	High	2	Yes	NA	BTT	HMII	dead	Died	5	26646	9754	36.61	98.69	58.32	7.57	64.88	35.12	2.06	0.55	46.50	12.62	17.70	23.18	15.07	53.31	12.50	83.40	1.10	15.64	1.76	28.12	0.00	7.54	16.46	0.55	1.92	6.04	13.00
96	22	67	older	Male	High	2	Yes	NA	BTT	HMII	dead	Died	8	73641	21944	29.80	96.79	48.91	4.47	62.74	37.26	2.84	3.16	37.47	17.79	24.95	19.79	19.46	38.11	9.46	79.79	1.89	13.58	1.24	27.05	20.00	14.11	13.37	0.21	2.74	10.53	3.31
97	23	68	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	0	44004	20612	46.84	43.65	15.61	21.54	49.64	50.36	10.37	0.57	15.94	45.15	33.18	5.73	24.69	5.68	20.25	44.58	1.55	36.17	0.88	67.23	27.27	32.35	15.48	20.02	4.49	10.94	2.76
98	23	68	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	1	99731	23424	23.49	91.81	41.33	10.89	35.24	64.76	11.83	1.62	18.71	59.21	16.19	5.89	28.16	5.89	18.53	39.98	4.66	42.16	2.33	84.07	2.63	22.55	15.21	40.28	6.07	9.70	3.96
99	23	68	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	8	129637	50001	38.57	90.65	44.41	6.05	30.66	69.34	17.12	4.05	10.62	66.50	19.60	3.28	26.41	3.46	13.37	39.74	7.15	39.49	11.91	71.97	2.70	32.48	33.14	43.91	7.81	22.70	3.82
100	23	68	older	Male	High	2	No	No	BTT	HMII	alive	Alive s/p OHT	14	134368	46363	34.50	88.32	34.25	4.63	32.47	67.53	12.86	8.38	8.91	64.89	22.88	3.32	24.78	3.89	13.58	36.16	6.01	40.33	13.77	72.11	1.89	29.89	26.36	42.80	6.17	15.97	3.29
101	24	64	older	Male	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	0	59881	8971	14.98	97.40	50.06	9.57	75.84	24.16	3.59	8.13	64.83	13.40	15.67	6.10	19.31	18.81	16.83	76.67	6.10	6.82	2.00	88.16	47.06	10.53	11.24	10.41	9.93	6.94	11.93
102	24	64	older	Male	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	1	60193	23874	39.66	96.72	58.83	16.24	80.37	19.63	1.87	3.31	68.58	10.62	14.72	6.08	21.33	22.83	3.80	80.61	3.63	8.03	0.33	94.16	57.26	2.96	5.01	10.19	4.72	1.41	1.24
103	24	64	older	Male	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	5	135335	27834	20.57	97.35	44.75	5.61	63.22	36.78	8.49	17.63	39.01	23.16	30.72	7.11	17.53	10.38	12.16	42.70	4.93	19.28	4.68	67.37	4.10	26.64	17.76	19.21	14.54	14.47	10.60
104	24	64	older	Male	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	8	247113	62042	25.11	96.41	43.93	14.82	78.42	21.58	2.52	5.24	65.35	12.22	17.53	4.90	20.39	17.20	2.93	77.86	2.11	8.78	2.94	84.64	9.70	10.64	13.93	9.90	11.59	8.48	3.51
105	24	64	older	Male	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	14	142053	56983	40.11	93.01	62.72	12.87	85.10	14.90	0.69	1.64	68.02	8.40	21.03	2.55	17.03	13.68	1.08	83.33	0.37	7.38	0.26	79.00	16.07	3.01	5.85	5.16	5.24	1.77	0.70
106	24	64	older	Male	High	2	No	Yes	BTT	HMII	alive	Alive s/p OHT	21	60384	17029	28.20	93.65	40.96	21.36	83.80	16.20	1.47	2.32	63.52	9.69	23.80	2.99	23.55	15.58	0.91	78.16	1.09	8.13	1.00	83.65	17.72	4.70	8.42	5.99	7.22	3.76	1.56
107	25	25	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	0	47645	37377	78.45	94.61	2.47	9.41	78.56	17.02	1.71	4.63	34.94	16.09	44.77	4.21	4.35	28.99	59.42	67.26	3.10	26.22	7.50	72.82	53.90	13.50	28.20	3.04	7.70	7.67	0.96
108	25	25	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	1	87452	43025	49.20	91.54	1.98	26.71	92.06	6.76	0.77	4.67	47.21	4.70	45.17	2.92	6.06	18.18	60.61	81.01	3.95	14.10	12.04	85.97	79.84	5.82	14.15	1.11	4.77	2.84	0.00
109	25	25	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	3	124056	80873	65.19	95.41	0.37	14.36	90.47	5.98	0.39	0.60	41.84	4.26	51.08	2.81	18.82	35.29	43.53	86.97	0.10	11.32	1.10	75.49	12.12	3.80	14.94	0.24	0.88	3.25	0.00
110	25	25	younger	Female	Low	3	No	Yes	BTT	HMII	alive	Alive s/p OHT	8	87840	57990	66.02	95.60	0.11	15.80	89.00	9.49	1.56	5.02	45.89	7.18	43.54	3.39	6.10	19.51	56.71	84.23	3.14	9.95	10.70	72.13	52.73	7.06	17.53	0.74	3.13	5.13	0.13
111	26	49	younger	Male	Low	3	No	NA	DT	HMII	alive	Alive	0	196315	122053	62.17	95.37	60.66	22.66	91.10	8.90	0.45	4.95	85.57	6.30	5.42	2.71	13.13	28.32	20.05	86.32	2.45	5.42	0.96	98.68	94.19	2.50	2.96	4.63	1.24	0.33	0.89
112	26	49	younger	Male	Low	3	No	NA	DT	HMII	alive	Alive	1	91909	37622	40.93	96.10	43.22	38.85	95.48	4.52	0.14	2.19	85.07	3.09	10.33	1.51	21.03	27.23	9.08	88.84	1.13	3.91	0.23	98.77	88.27	1.16	1.16	2.34	0.38	0.21	0.45
113	26	49	younger	Male	Low	3	No	NA	DT	HMII	alive	Alive	5	153762	86281	56.11	95.83	65.54	24.24	91.12	8.88	0.65	2.84	81.43	6.30	9.62	2.65	21.14	28.73	15.89	83.00	1.97	6.74	0.36	98.22	93.16	3.20	2.29	5.47	1.11	0.77	1.56
114	26	49	younger	Male	Low	3	No	NA	DT	HMII	alive	Alive	8	82406	47485	57.62	92.71	65.75	17.46	91.79	8.21	1.08	1.99	76.35	5.70	15.41	2.54	20.56	26.95	11.99	78.67	1.35	7.97	0.28	97.55	85.62	2.65	2.56	5.19	1.18	0.86	1.79
115	27	62	older	Male	High	2	No	NA	BTT	CMAG	dead	Died	0	105070	38521	36.66	98.39	47.84	12.88	60.82	39.18	11.08	2.89	33.87	29.43	26.68	10.01	19.72	21.92	2.99	59.35	14.64	20.23	3.45	81.40	46.10	24.47	24.27	24.04	11.00	17.02	1.51
116	27	62	older	Male	High	2	No	NA	BTT	CMAG	dead	Died	1	162336	73080	45.02	99.08	49.35	13.95	74.24	25.76	7.77	3.84	45.21	19.33	28.85	6.61	19.16	23.45	0.81	65.34	14.40	13.50	3.17	86.68	71.39	23.26	17.70	15.21	6.93	13.29	0.88
117	27	62	older	Male	High	2	No	NA	BTT	CMAG	dead	Died	5	111264	60487	54.36	99.15	45.94	16.13	73.16	26.84	4.64	1.82	52.04	20.07	20.94	6.95	22.93	23.47	0.50	77.11	9.30	9.58	1.20	88.92	81.82	22.20	19.01	16.98	6.89	13.57	1.36
118	27	62	older	Male	High	2	No	NA	BTT	CMAG	dead	Died	8	165870	67935	40.96	98.88	43.01	15.86	70.23	29.77	7.85	1.88	40.39	20.49	29.66	9.46	22.70	29.00	1.88	68.61	9.60	16.11	1.47	89.57	69.00	20.12	16.59	17.10	7.02	12.40	2.06
119	28	72	older	Male	Low	4	Yes	NA	BTT	HMII	dead	Died	3	424375	21971	5.18	71.78	56.98	18.53	79.84	20.16	3.08	0.51	57.91	14.07	24.06	3.97	31.74	15.95	13.65	57.94	0.65	28.51	0.17	93.57	53.33	2.87	9.69	15.61	3.87	2.16	3.99
120	28	72	older	Male	Low	4	Yes	NA	BTT	HMII	dead	Died	8	177249	103442	58.36	76.40	15.43	5.64	65.13	34.87	5.21	0.20	35.56	27.91	31.90	4.62	26.60	8.28	13.24	67.89	0.29	20.13	0.00	78.60	22.22	5.50	16.09	25.56	6.80	4.73	0.89

Identification of clusters

We double standardized the data and biclustered it. We found 3 clusters of patients, and 3 clusters of B-cells.

require(pheatmap, quietly = T)
require(RColorBrewer, quietly = T)


colorfun <- function(nlevels, ...){
    if(nlevels > 10) return(standardColors(nlevels))
    if(nlevels <=10) return(pal_d3(...)(nlevels))
} 

annotation.col <- df.raw[,c(1,3,4,6,7,8,9,10,11,12)]
annotation.colors <- lapply(colnames(annotation.col), function(nam){
    colrs <- colorfun(nlevels(annotation.col[[nam]]))
    names(colrs) <- levels(annotation.col[[nam]])
    return(colrs)
})
names(annotation.colors) <- colnames(annotation.col)

df.patient <- double_standardize(t(df.raw[,c(2,14:42)]))
pheatmap(df.patient, 
         color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdBu")))(100),
         breaks = c(seq(-max(abs(df.patient), na.rm=T), 0, length.out = 50), seq(0.001, max(abs(df.patient), na.rm=T), length.out = 50)),
         #scale = "row",
         fontsize = 7,
         cutree_rows = 3,
         cutree_cols = 3,
         annotation_col = annotation.col,
         annotation_colors = annotation.colors
)

Metadata associations

We computed the statistical associations between sample groups specified by the metadata.

Survival

We computed metadata factors that were statistically associated with survival.

df.factors <- df[match(levels(df$PatientID), df$PatientID),]
isfactor <- which(sapply(df.factors, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("Survival","Outcome", "LowIntermacs"))], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$Survival)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)

	OR	lowerCL	upperCL	pvalue	qvalue	Factor
1	4.900	0.7440	39.00	0.0974	0.681	RVAD
2	NA	NA	NA	0.354	0.823	Device Type
3	2.610	0.1070	188.00	0.569	0.823	Sensitized
4	NA	NA	NA	0.616	0.823	InterMACS
5	0.677	0.0853	5.94	0.674	0.823	Sex
6	1.480	0.2450	9.82	0.705	0.823	AgeGreater60
7	0.879	0.0656	7.89	1	1	VAD Indication

InterMACS

We computed metadata factors that were statistically associated with a binarized InterMACS score.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("LowIntermacs", "InterMACS"))], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$LowIntermacs)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)

	OR	lowerCL	upperCL	pvalue	qvalue	Factor
1	NA	NA	NA	0.135	0.983	Device Type
2	3.070	0.520	21.10	0.246	0.983	AgeGreater60
3	0.583	0.062	5.43	0.653	1	VAD Indication
4	NA	NA	NA	0.689	1	Outcome
5	1.830	0.288	14.60	0.689	1	RVAD
6	1.310	0.116	15.90	1	1	Sensitized
7	1.210	0.140	9.40	1	1	Sex
8	0.956	0.153	6.40	1	1	Survival

Sensitization

We computed metadata factors that were statistically associated with sensitization.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("Sensitized"))], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$Sensitized)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)

	OR	lowerCL	upperCL	pvalue	qvalue	Factor
1	0.127	0.00192	1.99	0.119	1	Sex
2	2.610	0.10700	188.00	0.569	1	Survival
3	NA	NA	NA	0.569	1	Outcome
4	0.426	0.02530	5.03	0.608	1	RVAD
5	NA	NA	NA	0.843	1	Device Type
6	1.310	0.11600	15.90	1	1	LowIntermacs
7	1.230	0.10100	15.40	1	1	AgeGreater60
8	NA	NA	NA	1	1	InterMACS
9	0.000	0.00000	Inf	1	1	VAD Indication

Sex

We computed metadata factors that were statistically associated with sex.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("Sex"))], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$Sex)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)

	OR	lowerCL	upperCL	pvalue	qvalue	Factor
1	0.127	0.00192	1.99	0.119	0.89	Sensitized
2	3.860	0.48600	49.80	0.198	0.89	AgeGreater60
3	NA	NA	NA	0.368	1	Outcome
4	NA	NA	NA	0.447	1	Device Type
5	0.677	0.08530	5.94	0.674	1	RVAD
6	0.677	0.08530	5.94	0.674	1	Survival
7	1.210	0.14000	9.40	1	1	LowIntermacs
8	NA	NA	NA	1	1	InterMACS
9	1.840	0.15200	103.00	1	1	VAD Indication

Age

We computed metadata factors that were statistically associated with age.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("AgeGreater60"))], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$AgeGreater60)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)

	OR	lowerCL	upperCL	pvalue	qvalue	Factor
1	0.228	0.0274	1.45	0.114	0.442	RVAD
2	NA	NA	NA	0.128	0.442	InterMACS
3	NA	NA	NA	0.159	0.442	Outcome
4	3.860	0.4860	49.80	0.198	0.442	Sex
5	3.070	0.5200	21.10	0.246	0.442	LowIntermacs
6	NA	NA	NA	0.35	0.525	Device Type
7	1.950	0.2230	25.90	0.655	0.794	VAD Indication
8	1.480	0.2450	9.82	0.705	0.794	Survival
9	1.230	0.1010	15.40	1	1	Sensitized

VAD Indication

We computed metadata factors that were statistically associated with VAD Indication.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("VAD Indication"))], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$`VAD Indication`)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)

	OR	lowerCL	upperCL	pvalue	qvalue	Factor
1	NA	NA	NA	0.000242	0.00217	Outcome
2	0.000	0.0000	1.34	0.0619	0.278	RVAD
3	0.583	0.0620	5.43	0.653	1	LowIntermacs
4	1.950	0.2230	25.90	0.655	1	AgeGreater60
5	NA	NA	NA	0.837	1	InterMACS
6	NA	NA	NA	0.877	1	Device Type
7	1.840	0.1520	103.00	1	1	Sex
8	0.000	0.0000	Inf	1	1	Sensitized
9	0.879	0.0656	7.89	1	1	Survival

RVAD

We computed metadata factors that were statistically associated with RVAD.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("RVAD"))], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$RVAD)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)

	OR	lowerCL	upperCL	pvalue	qvalue	Factor
1	NA	NA	NA	0.00868	0.0782	Device Type
2	0.000	0.0000	1.34	0.0619	0.205	VAD Indication
3	NA	NA	NA	0.0739	0.205	Outcome
4	4.900	0.7440	39.00	0.0974	0.205	Survival
5	0.228	0.0274	1.45	0.114	0.205	AgeGreater60
6	NA	NA	NA	0.327	0.49	InterMACS
7	0.426	0.0253	5.03	0.608	0.689	Sensitized
8	0.677	0.0853	5.94	0.674	0.689	Sex
9	1.830	0.2880	14.60	0.689	0.689	LowIntermacs

Device-type

We computed metadata factors that were statistically associated with Device Type.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("Device Type"))], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$`Device Type`)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = T, B = 10000)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = T, B = 10000)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)

	OR	lowerCL	upperCL	pvalue	qvalue	Factor
1	NA	NA	NA	0.0089	0.0801	RVAD
2	NA	NA	NA	0.139	0.528	LowIntermacs
3	NA	NA	NA	0.274	0.528	InterMACS
4	NA	NA	NA	0.347	0.528	Outcome
5	NA	NA	NA	0.352	0.528	AgeGreater60
6	NA	NA	NA	0.352	0.528	Survival
7	NA	NA	NA	0.443	0.569	Sex
8	NA	NA	NA	0.835	0.874	Sensitized
9	NA	NA	NA	0.874	0.874	VAD Indication

Visualization of metadata associations

We used the GGally package to visualize pairwise dependencies between sample groups, separating out survivors and non-survivors (indicated by color).

df.temp <- df.factors
colnames(df.temp) <- make.names(colnames(df), unique = T)
#invisible(
suppressWarnings(
    suppressMessages(
        ggpairs(df.temp, 
                #mapping = aes(color = Survival), 
                columns = colnames(df.temp)[c(2,3,4,5,6,7,8,9,10,11,12)]) 
    ))

#)

Variability of MCS devices

PCA

Using Principal Component Analysis (PCA), we saw large variability between device types, compared to the variability within device types. We also saw large variability between individual patients. None of the other features were clearly separable.

suppressMessages(require(ggbiplot, quietly = T))
require(ggsci, quietly = T)

isna <- unique(unlist(apply(df[,14:42], 2, function(x) which(is.na(x)))))
pca <- prcomp(double_standardize(df[-isna, 14:42]), center = TRUE, scale. = TRUE)

colorfun <- function(grouping, ...){
    if(nlevels(factor(grouping)) > 10) scale_color_discrete(...)
    if(nlevels(factor(grouping)) <=10) scale_color_d3(...)
} 

plotvars <- names(df)[c(10,1:9,11:17)]
plots.pca <- mclapply(plotvars, function(this_var){
    this_groups <- df[-isna, this_var]
    if(is.numeric(this_groups)) this_groups <- cut(this_groups, 4)
    
    ggbiplot(pca,
             groups = this_groups,
             ellipse = TRUE,
             alpha = 0.3,
             varname.size = 1.2) + 
        colorfun(this_groups, name = this_var) +
        ggtitle(paste0(this_var, " variability")) +
        theme_classic() 
}, mc.cores = detectCores()-1)
names(plots.pca) <- plotvars
#plots.pca$PatientID
#plots.pca$`Device Type`

for(ii in 1:length(plots.pca)){
    cat("  \n###", names(plots.pca)[ii], "\n")
    suppressWarnings(print(plots.pca[[ii]]))
    cat("  \n")
}

Device Type

PatientID

Age

AgeGreater60

Sex

LowIntermacs

InterMACS

RVAD

Sensitized

VAD Indication

Survival

Outcome

Time

num Total PBMC

num lymph

lymph

live lymph

Fisher’s exact test

We analyzed the dependency of device types on the other discrete variables using Fisher’s exact test.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[names(isfactor) != "Device Type"], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$`Device Type`)
                      })

fisher.p <- sapply(contingency, function(this_table){
    fisher.test(this_table, simulate.p.value = T, B = 10000)$p
})
fisher.q <- p.adjust(fisher.p, method = "BH")

signif.ix <- which(fisher.p < 0.05)
signif.order <- sort(fisher.q[signif.ix], index.return = T)$ix
for(this_ix in signif.ix[signif.order]){
    cat("  \n###", names(contingency)[this_ix], "\n")
    cat(paste0("Benjamini-Hochberg qvalue = ", signif(fisher.q[this_ix], 2)),
        ". \n")
    print(kable(contingency[this_ix][[1]],  row.names = T) %>%
              kable_styling(bootstrap_options = c("striped", 
                                                  "hover", 
                                                  "condensed",
                                                  "responsive"),
                            font_size = 12)
    )
    cat("  \n")
}

RVAD

Benjamini-Hochberg qvalue = 0.077 .

	HMII	CMAG	HVAD	PVAD	TAH
No	15	1	2	0	0
Yes	5	0	0	3	2

One-way repeated measures anova

We first analyzed the differences in B-cell levels across device types using a one-way repeated measures anova. Here we report any variables that had a statistically significant variance (\(p<0.05\)) across devices, time, or their interaction.

suppressMessages(require(lmerTest, quietly = T))
suppressMessages(require(car, quietly = TRUE))

df.device <- df
colnames(df.device) <- make.names(colnames(df), unique = T)

varnames <- colnames(df.device)[14:42]
models.device <- mclapply(varnames, function(this_var){
    this_formula <- as.formula(paste0(this_var, " ~ Device.Type + (1|PatientID)"))
    invisible(suppressMessages(this_model <- lmer(this_formula, data = df.device)))
    this_anova <- Anova(this_model, type = 2)
    pvals <- this_anova$`Pr(>Chisq)`[c(1)]
    return(list(model = this_model,
                pvals = pvals))
}, mc.cores = detectCores()-1)
names(models.device) <- colnames(df)[14:42]

pvals <- do.call(rbind, lapply(models.device, function(x) x$pvals))
rownames(pvals) <- colnames(df)[14:42]
colnames(pvals) <- c("pvalue")

qBH <- matrix(p.adjust(pvals, method = "BH"), nrow = nrow(pvals))
rownames(qBH) <- colnames(df)[14:42]
colnames(qBH) <- c("qvalue")

sigvars <- apply(apply(pvals, 2, function(x) x<=0.05), 1, function(x) sum(x) > 0)

sig.qtable <- cbind(pvals,qBH)[sigvars,,drop=F][order(apply(pvals[sigvars,,drop=F], 1, min)),,drop=F]

qtable <- as.data.frame(signif(sig.qtable, 2))
qtable %>%
    mutate(
        `B-cell` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>%
    kable( escape = F,
           digits = 20,
           row.names = T,
           caption = "Significant device-type q-values") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  font_size = 12) %>%
    scroll_box(width = "100%")

Significant device-type q-values

	pvalue	qvalue	B-cell
1	0.011	0.31	CD27-IgD+ mature naive
2	0.045	0.45	CD19+CD5+

We plotted mean levels across time for each of the B-cells that showed a statistically significant effect across devices in the above mixed effect models. We drew attention to specific features that induced the positive test result, by listing the model parameters with \(p<0.05\) in the fit. Note that the reference level for the devices is the HeartMate-II (HMII).

require(reshape2, quietly = T)
df.long <- melt(df, id.vars = colnames(df)[1:13])
groups <- make.names(c("Device Type"))
names(df.long) <- make.names(names(df.long))

invisible(suppressMessages(require(Hmisc, quietly = T)))
stat_sum_df <- function(fun, geom="errorbar", ...) {
    stat_summary(fun.data = fun, geom = geom, width = 1, ...)
}

plots.ts <- mclapply(rownames(qtable), function(this_var){
    lapply(groups, function(this_groups){
        ggplot(subset(df.long, df.long$variable == this_var)) +
            aes(x = Time, y = value, group = PatientID) +
            aes_string(color = this_groups, fill = this_groups) +
            geom_line(alpha = 0) + 
            geom_point(alpha = 0) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("point"), position = position_dodge(.5)) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("line"), position = position_dodge(.5)) + 
            stat_sum_df(function(x) mean_cl_normal(x, conf.int = 0.68), mapping = aes_string(group = this_groups), position = position_dodge(.5)) + 
            #stat_smooth(aes_string(group = this_groups), method = "loess", span = 1) +
            scale_color_d3() + scale_fill_d3() +
            xlab("Time (days after surgery)") +
            ylab(this_var) + 
            ggtitle(paste(this_var)) +
            theme_classic()
    })
}, mc.cores = detectCores()-1)
names(plots.ts) <- rownames(qtable)

for(ii in c(1:length(plots.ts))){
    for(jj in 1:length(plots.ts[[ii]])){
        sumtable <- suppressMessages(summary(models.device[[rownames(qtable)[ii]]]$model))
        sumtable <- as.data.frame(sumtable$coefficients)[-1, ,drop=F] # drop intercept
        if(!("Pr(>|t|)" %in% colnames(sumtable))) next()
        sigsum <- sumtable[sumtable[,"Pr(>|t|)"] <= 0.05, , drop = F]
        cat("  \n###", rownames(qtable)[ii], "\n")
        print(kable(sigsum[order(sigsum[,"Pr(>|t|)"]),,drop=F], row.names = T) %>%
                  kable_styling(bootstrap_options = c("striped", 
                                                      "hover", 
                                                      "condensed",
                                                      "responsive"),
                                font_size = 12)
        )
        cat("  \n")
        suppressWarnings(print(plots.ts[[ii]][[jj]]))
        cat("  \n")
    }
}

CD27-IgD+ mature naive

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypePVAD	-33.85486	10.06687	23.54824	-3.362998	0.0026287

CD19+CD5+

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypePVAD	23.16366	8.615765	23.18582	2.688521	0.0130627

Linar mixed-effect model

We next analyzed the differences in B-cell levels across device types using a linear mixed effect model with time as a continuous variable, and included the interaction term. Here we report variables that had a statistically significant variance (Benjamini-Hochberg \(p<0.05\)) across devices, time, or their interaction.

suppressMessages(require(lmerTest, quietly = T))
suppressMessages(require(car, quietly = TRUE))

df.device <- df
colnames(df.device) <- make.names(colnames(df), unique = T)

varnames <- colnames(df.device)[14:42]
models.device <- mclapply(varnames, function(this_var){
    this_formula <- as.formula(paste0(this_var, " ~ Device.Type * Time + (1|PatientID)"))
    invisible(suppressMessages(this_model <- lmer(this_formula, data = df.device)))
    this_anova <- Anova(this_model, type = 2)
    pvals <- this_anova$`Pr(>Chisq)`[c(1,2,3)]
    return(list(model = this_model,
                pvals = pvals))
}, mc.cores = detectCores()-1)
names(models.device) <- colnames(df)[14:42]

pvals <- do.call(rbind, lapply(models.device, function(x) x$pvals))
rownames(pvals) <- varnames
colnames(pvals) <- c("Device", "Time", "Device:Time")

qBH <- matrix(p.adjust(pvals, method = "BH"), nrow = nrow(pvals))
rownames(qBH) <- colnames(df)[14:42]
colnames(qBH) <- colnames(pvals)

sigvars <- apply(apply(pvals, 2, function(x) x<=0.05), 1, function(x) sum(x) > 0)
sig.qtable <- qBH[sigvars,][order(apply(qBH[sigvars,], 1, min)),]

qtable <- as.data.frame(signif(sig.qtable, 2))
qtable %>%
    mutate(
        `B-cell` = row.names(.),
        Device = cell_spec(Device, color = ifelse(qtable$Device > 0.05, "grey", "red")),
        `Device:Time` = cell_spec(`Device:Time`, color = ifelse(qtable$`Device:Time` > 0.05, "grey", "red")),
        `Time` = cell_spec(`Time`, color = ifelse(qtable$`Time` > 0.05, "grey", "red"))
    ) %>%
    kable( escape = F,
           digits = 20,
           row.names = T,
           caption = "Significant device-type q-values") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  font_size = 12) %>%
    scroll_box(width = "100%")

Significant device-type q-values

	Device	Time	Device:Time	B-cell
1	0.6	0.0023	0.007	CD19+CD11b+
2	0.5	0.025	0.82	CD268 of +27-38++transitional
3	0.82	0.025	0.82	CD19+CD5+CD11b+
4	0.85	0.089	0.82	CD19+27-38+CD5+transitionals
5	0.4	0.089	0.87	CD19+CD268+
6	0.85	0.093	0.96	CD19 of live lymph
7	0.47	0.24	0.093	CD27+38++plasma blasts
8	0.15	0.2	0.24	CD27-IgD+ mature naive
9	0.47	0.18	0.82	CD27+IgD- switched memory
10	0.26	0.19	0.93	CD19+CD5+
11	0.82	0.19	0.87	CD19+CD5+CD24hi
12	0.82	0.2	0.82	num lymph

We plotted mean levels across time for each of the B-cells that showed a statistically significant effect across devices in the above mixed effect models. We drew attention to specific features that induced the positive test result, by listing the model parameters with \(p<0.05\) in the multivariate fit. Note that the reference level for the devices is the HeartMate-II (HMII).

require(reshape2, quietly = T)
df.long <- melt(df, id.vars = colnames(df)[1:13])
groups <- make.names(c("Device Type"))
names(df.long) <- make.names(names(df.long))

invisible(suppressMessages(require(Hmisc, quietly = T)))
stat_sum_df <- function(fun, geom="errorbar", ...) {
    stat_summary(fun.data = fun, geom = geom, width = 1, ...)
}

plots.ts <- mclapply(rownames(qtable), function(this_var){
    lapply(groups, function(this_groups){
        ggplot(subset(df.long, df.long$variable == this_var)) +
            aes(x = Time, y = value, group = PatientID) +
            aes_string(color = this_groups, fill = this_groups) +
            geom_line(alpha = 0) + 
            geom_point(alpha = 0) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("point"), position = position_dodge(.5)) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("line"), position = position_dodge(.5)) + 
            stat_sum_df(function(x) mean_cl_normal(x, conf.int = 0.68), mapping = aes_string(group = this_groups), position = position_dodge(.5)) + 
            #stat_smooth(aes_string(group = this_groups), method = "loess", span = 1) +
            scale_color_d3() + scale_fill_d3() +
            xlab("Time (days after surgery)") +
            ylab(this_var) + 
            ggtitle(paste(this_var)) +
            theme_classic()
    })
}, mc.cores = detectCores()-1)
names(plots.ts) <- rownames(qtable)

for(ii in c(1:length(plots.ts))){
    for(jj in 1:length(plots.ts[[ii]])){
        sumtable <- suppressMessages(summary(models.device[[rownames(qtable)[ii]]]$model))
        sumtable <- as.data.frame(sumtable$coefficients)[-1, ,drop=F] # drop intercept
        if(!("Pr(>|t|)" %in% colnames(sumtable))) next()
        sigsum <- sumtable[sumtable[,"Pr(>|t|)"] <= 0.05, , drop = F]
        cat("  \n###", rownames(qtable)[ii], "\n")
        print(kable(sigsum[order(sigsum[,"Pr(>|t|)"]),,drop=F], row.names = T) %>%
                  kable_styling(bootstrap_options = c("striped", 
                                                      "hover", 
                                                      "condensed",
                                                      "responsive"),
                                font_size = 12)
        )
        cat("  \n")
        suppressWarnings(print(plots.ts[[ii]][[jj]]))
        cat("  \n")
    }
}

CD19+CD11b+

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypeTAH:Time	2.1101624	0.5538816	89.91928	3.809771	0.0002541
Device.TypePVAD:Time	0.8054722	0.2943289	88.29844	2.736640	0.0075041
Time	0.2628129	0.1188074	89.80591	2.212092	0.0294973

CD268 of +27-38++transitional

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	-0.7873216	0.3030382	89.7994	-2.598094	0.0109556

CD19+CD5+CD11b+

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	0.290629	0.1150892	91.27507	2.525249	0.0132869

CD19+27-38+CD5+transitionals

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	0.0910247	0.0395736	89.94153	2.30014	0.0237538

CD19+CD268+

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	-0.6425061	0.230173	88.33500	-2.791405	0.0064310
Device.TypePVAD	-37.7748686	16.138544	27.53028	-2.340662	0.0267323

CD19 of live lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	-0.3683497	0.140373	93.29159	-2.624078	0.0101509

CD27+38++plasma blasts

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypePVAD:Time	0.2768179	0.1033923	88.20169	2.677355	0.0088487

CD27-IgD+ mature naive

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypePVAD	-32.958219	10.7368535	29.72852	-3.069635	0.0045463
Device.TypeTAH:Time	-2.247904	0.8166692	89.34072	-2.752527	0.0071618

CD27+IgD- switched memory

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	0.3654018	0.1606941	88.62502	2.273897	0.0253875
Device.TypePVAD	25.0335542	10.8348873	28.15496	2.310458	0.0284029

CD19+CD5+

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypePVAD	19.97858	9.433345	33.29567	2.117868	0.0417364

CD19+CD5+CD24hi

	Estimate	Std. Error	df	t value	Pr(>|t|)

Two-way repeated measures anova

Finally, we attempted to analyze the differences in B-cell levels across device types using a two-way repeated measures ANOVA. Here we report variables that had a statistically significant variance (\(p<0.05\)) across devices, times, or their interaction. As there were 5 device types and 7 timepoints, but only 120 samples, this model is severely underpowered. The posterior belief in any of these results should therefore be quite small (as a consequence of Bayes rule).

suppressMessages(require(lmerTest, quietly = T))
suppressMessages(require(car, quietly = TRUE))

df.device <- df
colnames(df.device) <- make.names(colnames(df), unique = T)

varnames <- colnames(df.device)[14:42]
models.device <- mclapply(varnames, function(this_var){
    this_formula <- as.formula(paste0(this_var, " ~ Device.Type * factor(Time) + (1|PatientID)"))
    invisible(suppressMessages(this_model <- lmer(this_formula, data = df.device)))
    this_anova <- Anova(this_model, type = 2)
    pvals <- this_anova$`Pr(>Chisq)`[c(1,2,3)]
    return(list(model = this_model,
                pvals = pvals))
}, mc.cores = detectCores()-1)
names(models.device) <- colnames(df)[14:42]

pvals <- do.call(rbind, lapply(models.device, function(x) x$pvals))
rownames(pvals) <- varnames
colnames(pvals) <- c("Device", "factor(Time)", "Device:factor(Time)")

qBH <- matrix(p.adjust(pvals, method = "BH"), nrow = nrow(pvals))
rownames(qBH) <- colnames(df)[14:42]
colnames(qBH) <- colnames(pvals)

sigvars <- apply(apply(pvals, 2, function(x) x<=0.05), 1, function(x) sum(x) > 0)
sig.qtable <- qBH[sigvars,][order(apply(qBH[sigvars,], 1, min)),]

qtable <- as.data.frame(signif(sig.qtable, 2))
qtable %>%
    mutate(
        `B-cell` = row.names(.),
        Device = cell_spec(Device, color = ifelse(qtable$Device > 0.05, "grey", "red")),
        `Device:factor(Time)` = cell_spec(`Device:factor(Time)`, color = ifelse(qtable$`Device:factor(Time)` > 0.05, "grey", "red")),
        `factor(Time)` = cell_spec(`factor(Time)`, color = ifelse(qtable$`factor(Time)` > 0.05, "grey", "red"))
    ) %>%
    kable( escape = F,
           digits = 20,
           row.names = T,
           caption = "Significant device-type q-values") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  font_size = 12) %>%
    scroll_box(width = "100%")

Significant device-type q-values

	Device	factor(Time)	Device:factor(Time)	B-cell
1	0.16	0.071	3e-10	CD19+CD5+
2	0.33	0.041	8.7e-06	CD27+IgD- switched memory
3	0.54	0.071	7.8e-05	CD19+CD27+
4	0.54	0.31	7.8e-05	CD27-IgD- switched memory
5	0.59	4e-04	0.0017	CD19+CD11b+
6	0.071	0.16	0.0034	CD27-IgD+ mature naive
7	0.93	0.0058	0.98	CD19 of live lymph
8	0.93	0.94	0.012	CD19CD24hiCD38-memory
9	0.27	0.041	1	CD27+IgD-IgM+ switched memory
10	0.98	0.93	0.057	CD19+24dim38dim naive mature
11	0.93	0.54	0.069	CD19+CD27+CD24hi
12	0.54	0.069	0.93	CD268 of +27-38++transitional
13	0.85	0.1	0.59	CD19+CD5+CD11b+
14	0.31	0.16	0.11	CD19+CD268+
15	0.44	0.22	0.16	CD27+38++plasma blasts

We plotted mean levels across time for each of the B-cells that showed a statistically significant effect across devices in the above mixed effect models. We drew attention to specific features that induced the positive test result, by listing the model parameters with \(p<0.05\) in the multivariate fit. Note that the reference level for the time comparisons is timepoint 0, and the reference level for the devices is the HeartMate-II (HMII).

require(reshape2, quietly = T)
df.long <- melt(df, id.vars = colnames(df)[1:13])
groups <- make.names(c("Device Type"))
names(df.long) <- make.names(names(df.long))

invisible(suppressMessages(require(Hmisc, quietly = T)))
stat_sum_df <- function(fun, geom="errorbar", ...) {
    stat_summary(fun.data = fun, geom = geom, width = 1, ...)
}

plots.ts <- mclapply(rownames(qtable), function(this_var){
    lapply(groups, function(this_groups){
        ggplot(subset(df.long, df.long$variable == this_var)) +
            aes(x = Time, y = value, group = PatientID) +
            aes_string(color = this_groups, fill = this_groups) +
            geom_line(alpha = 0) + 
            geom_point(alpha = 0) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("point"), position = position_dodge(.5)) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("line"), position = position_dodge(.5)) + 
            stat_sum_df(function(x) mean_cl_normal(x, conf.int = 0.68), mapping = aes_string(group = this_groups), position = position_dodge(.5)) + 
            #stat_smooth(aes_string(group = this_groups), method = "loess", span = 1) +
            scale_color_d3() + scale_fill_d3() +
            xlab("Time (days after surgery)") +
            ylab(this_var) + 
            ggtitle(paste(this_var)) +
            theme_classic()
    })
}, mc.cores = detectCores()-1)
names(plots.ts) <- rownames(qtable)

for(ii in c(1:length(plots.ts))){
    for(jj in 1:length(plots.ts[[ii]])){
        sumtable <- suppressMessages(summary(models.device[[rownames(qtable)[ii]]]$model))
        sumtable <- as.data.frame(sumtable$coefficients)[-1, ,drop=F] # drop intercept
        if(!("Pr(>|t|)" %in% colnames(sumtable))) next()
        sigsum <- sumtable[sumtable[,"Pr(>|t|)"] <= 0.05, , drop = F]
        cat("  \n###", rownames(qtable)[ii], "\n")
        print(kable(sigsum[order(sigsum[,"Pr(>|t|)"]),], row.names = T) %>%
                  kable_styling(bootstrap_options = c("striped", 
                                                      "hover", 
                                                      "condensed",
                                                      "responsive"),
                                font_size = 12)
        )
        cat("  \n")
        suppressWarnings(print(plots.ts[[ii]][[jj]]))
        cat("  \n")
    }
}

CD19+CD5+

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypeTAH:factor(Time)5	79.25674	11.074689	67.92369	7.156565	0.0000000
Device.TypePVAD	35.52841	9.894557	43.53668	3.590703	0.0008326
Device.TypePVAD:factor(Time)3	-25.18076	8.751886	71.43140	-2.877181	0.0052873
Device.TypePVAD:factor(Time)5	-31.58267	11.241069	73.06105	-2.809579	0.0063611
Device.TypePVAD:factor(Time)8	-20.60508	7.572762	69.16187	-2.720946	0.0082277

CD27+IgD- switched memory

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypePVAD:factor(Time)3	-44.11347	8.159997	70.32575	-5.406064	0.0000008
Device.TypePVAD	46.55839	10.997051	35.92281	4.233715	0.0001521
Device.TypePVAD:factor(Time)8	-25.13824	7.037695	68.87488	-3.571943	0.0006518
Device.TypePVAD:factor(Time)5	-35.76427	10.505189	71.34332	-3.404439	0.0010914
Device.TypeHVAD:factor(Time)1	22.95446	9.340788	68.84101	2.457443	0.0165160
Device.TypeTAH:factor(Time)1	19.80313	9.340788	68.84101	2.120071	0.0376068

CD19+CD27+

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypePVAD:factor(Time)3	-44.86147	9.352169	70.81469	-4.796905	0.0000087
Device.TypePVAD	40.90474	11.677321	38.82945	3.502922	0.0011749
Device.TypePVAD:factor(Time)5	-36.38395	12.029880	72.05165	-3.024465	0.0034502
Device.TypePVAD:factor(Time)8	-22.52547	8.075489	69.06610	-2.789363	0.0068193
Device.TypeHVAD:factor(Time)1	26.28504	10.718455	69.03694	2.452316	0.0167262

CD27-IgD- switched memory

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypeTAH:factor(Time)5	71.68679	13.81377	66.76087	5.189516	0.0000021
Device.TypeTAH:factor(Time)14	34.72518	15.34586	76.65452	2.262837	0.0264767
Device.TypePVAD:factor(Time)3	-22.85502	10.65984	75.69337	-2.144030	0.0352412

CD19+CD11b+

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypeTAH:factor(Time)14	52.95791	10.085405	72.02616	5.250946	0.0000015
Device.TypeTAH:factor(Time)1	21.89638	8.043661	70.05677	2.722190	0.0081768
Device.TypePVAD:factor(Time)14	17.80531	7.197456	72.09536	2.473834	0.0157228
Device.TypeTAH:factor(Time)3	16.76692	8.048562	70.13474	2.083219	0.0408780

CD27-IgD+ mature naive

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypeTAH:factor(Time)5	-49.338009	12.632153	68.31609	-3.905748	0.0002179
Device.TypePVAD	-40.588774	11.742845	41.95252	-3.456469	0.0012664
Device.TypeTAH:factor(Time)14	-42.518754	14.393799	71.18856	-2.953963	0.0042493
Device.TypeTAH:factor(Time)1	-29.588950	11.470396	69.40423	-2.579593	0.0120122
factor(Time)8	-7.421838	2.996063	68.79529	-2.477197	0.0157031
Device.TypeHVAD:factor(Time)1	-23.895494	11.470396	69.40423	-2.083232	0.0409155

CD19 of live lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)
factor(Time)1	8.670607	2.838543	69.69728	3.054598	0.0031913
Device.TypeHVAD:factor(Time)5	18.621453	8.106642	68.64256	2.297061	0.0246707
Device.TypeHVAD:factor(Time)3	16.313101	7.998551	68.60746	2.039507	0.0452529

CD19CD24hiCD38-memory

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypeTAH:factor(Time)5	-22.909686	6.066396	68.15378	-3.776490	0.0003363
Device.TypeHVAD:factor(Time)1	12.942595	5.519995	68.69011	2.344675	0.0219407
factor(Time)1	-3.388861	1.568014	68.34505	-2.161244	0.0341830

CD27+IgD-IgM+ switched memory

	Estimate	Std. Error	df	t value	Pr(>|t|)
factor(Time)5	10.711338	3.308708	71.05832	3.237317	0.0018346
factor(Time)3	8.467757	3.003200	71.42351	2.819578	0.0062189
factor(Time)14	8.273763	3.705650	70.67012	2.232742	0.0287354

CD19+24dim38dim naive mature

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypePVAD:factor(Time)3	-51.79949	13.95422	76.12859	-3.712101	0.0003892

CD19+CD27+CD24hi

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypeTAH:factor(Time)5	-22.54139	9.215036	68.32836	-2.446154	0.0170182
factor(Time)8	4.88642	2.184243	68.95972	2.237123	0.0285135
Device.TypeHVAD:factor(Time)1	18.31135	8.355636	69.79036	2.191497	0.0317555

CD268 of +27-38++transitional

	Estimate	Std. Error	df	t value	Pr(>|t|)
factor(Time)8	-13.21005	5.828370	69.05529	-2.266508	0.0265563
factor(Time)14	-16.08481	7.874549	69.65369	-2.042632	0.0448754

CD19+CD5+CD11b+

	Estimate	Std. Error	df	t value	Pr(>|t|)

CD19+CD268+

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypeTAH:factor(Time)5	-66.88010	17.16044	67.82057	-3.897342	0.0002253
Device.TypePVAD	-39.35919	17.41639	37.46865	-2.259893	0.0297274

CD27+38++plasma blasts

	Estimate	Std. Error	df	t value	Pr(>|t|)
Device.TypeTAH:factor(Time)1	7.049114	2.979632	70.86462	2.365767	0.0207319
Device.TypeCMAG	7.332777	3.630995	44.46774	2.019495	0.0494876

HeartMate-II analysis

The HeartMate-II (HMII) recipients were the largest group, and we analyzed them by themselves due to the previously observed variability across devices.

require(reshape2, quietly = T)
df.HMII <- subset(df, df$`Device Type`=="HMII")
df.long <- melt(df.HMII, id.vars = colnames(df)[1:13])

groups <- make.names(c("AgeGreater60", 
                       "Sex",
                       "LowIntermacs",
                       "RVAD", 
                       "Sensitized",
                       "VAD Indication", 
                       #"Device Type", 
                       "Survival",
                       "Outcome"))

names(df.long) <- make.names(names(df.long))

plots.ts <- mclapply(unique(df.long$variable), function(this_var){
    lapply(groups, function(this_groups){
        ggplot(subset(df.long, df.long$variable == this_var)) +
            aes(x = Time, y = value, group = PatientID) +
            aes_string(color = this_groups, fill = this_groups) +
            geom_line(alpha = 0) + 
            geom_point(alpha = 0) + 
            #stat_smooth(aes_string(group = this_groups), method = "loess", span = 1) +
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("point"), position = position_dodge(.5)) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("line"), position = position_dodge(.5)) + 
            stat_sum_df(function(x) mean_cl_normal(x, conf.int = 0.68), mapping = aes_string(group = this_groups), position = position_dodge(.5)) + 
            scale_color_aaas() + scale_fill_aaas() +
            xlab("Time (days after surgery)") +
            ylab(this_var) + 
            ggtitle(paste(this_var)) +
            theme_classic()
    })
}, mc.cores = detectCores()-1)
names(plots.ts) <- unique(df.long$variable)

# for(ii in c(1:length(plots.ts))){
#     for(jj in 1:length(plots.ts[[ii]])){
#         suppressWarnings(print(plots.ts[[ii]][[jj]]))
#     }
# } 

PCA

Using PCA, we found large variability between individual patients, compared to the variability within individual patients. None of the other features were clearly separable.

suppressMessages(require(ggbiplot, quietly = T))
require(ggsci, quietly = T)

# Efron's double standardization
double_standardize <- function(x, niter = 100) {
    for(i in 1:niter) x <- t(scale(t(scale(x))))
    return(as.data.frame(x))
}

isna <- unique(unlist(apply(df.HMII[,14:42], 2, function(x) which(is.na(x)))))
pca <- prcomp(double_standardize(df.HMII[-isna, 14:42]), center = TRUE, scale. = TRUE)

colorfun <- function(grouping, ...){
    if(nlevels(factor(grouping)) > 10) scale_color_discrete(...)
    if(nlevels(factor(grouping)) <=10) scale_color_d3(...)
} 

plots.pca <- mclapply(names(df.HMII)[1:17], function(this_var){
    this_groups <- df.HMII[-isna, this_var]
    if(is.numeric(this_groups)) this_groups <- cut(this_groups, 4)
    
    ggbiplot(pca,
             groups = this_groups,
             ellipse = TRUE,
             alpha = 0.3,
             varname.size = 1.2) + 
        colorfun(this_groups, name = this_var) +
        ggtitle(paste0(this_var, " variability")) +
        theme_classic() 
}, mc.cores = detectCores()-1)

names(plots.pca) <- names(df.HMII)[1:17]
#plots.pca$PatientID
#plots.pca$`Device Type`

for(ii in 1:length(plots.pca)){
    cat("  \n###", names(plots.pca)[ii], "\n")
    suppressWarnings(print(plots.pca[[ii]]))
    cat("  \n")
}

PatientID

Age

AgeGreater60

Sex

LowIntermacs

InterMACS

RVAD

Sensitized

VAD Indication

Device Type

Survival

Outcome

Time

num Total PBMC

num lymph

lymph

live lymph

One-way repeated measures anova

We analyzed the differences in B-cell levels for various features using a one-way repeated measures anova. Here we report variables that had a statistically significant variance (\(p<0.05\)) across groups, or groups at each timepoint.

suppressMessages(require(lmerTest, quietly = TRUE))
suppressMessages(require(car, quietly = TRUE))
require(reshape2, quietly = TRUE)

df.lmer <- df.HMII
names(df.lmer) <- make.names(names(df.lmer), unique = TRUE)

groupvars.ix <- c(3,4,5,7,8,9,11)
groupvars <- names(df.lmer)[groupvars.ix]

bcells.ix <- c(14:42)
bcells <- names(df.lmer)[bcells.ix]

models.b <- mclapply(groupvars, function(this_groupvar){
    models.bcells <- lapply(bcells, function(this_bcell){
        this_formula <- as.formula(paste0(this_bcell, " ~ ", this_groupvar, 
                                          " + (1|PatientID)"))
        suppressMessages(suppressWarnings(this_model <- lmer(this_formula, data = droplevels(df.lmer))))
        this_anova <- Anova(this_model, type = 2)
        this_pvalues <- this_anova$`Pr(>Chisq)`
        names(this_pvalues) <- rownames(this_anova)
        #return(this_pvalues)
        return(list(model = this_model,
                    pvals = this_pvalues))
    })
    names(models.bcells) <- colnames(df)[14:42]
    pvalues <- do.call(rbind, lapply(models.bcells, function(x) x$pvals))
    rownames(pvalues) <- bcells
    #return(pvalues)
    return(list(model = models.bcells,
                pvals = pvalues))
}, mc.cores = detectCores()-1)
names(models.b) <- groupvars
pvals <- lapply(models.b, function(x) x$pvals)
# something wrong here
names(pvals) <- groupvars
pvals.matrix <- do.call(cbind, lapply(pvals, function(this_pval) this_pval[,c(1)]))



# Benjamini Hochberg
# qBH <- matrix(p.adjust(as.numeric(pvals.matrix), 
#                        method = "BH"), 
#               nrow = nrow(pvals.matrix), 
#               ncol = ncol(pvals.matrix), 
#               byrow = F) 
# rownames(qBH) <- rownames(pvals.matrix)
# colnames(qBH) <- colnames(pvals.matrix)
# rownames(qBH) <- names(df)[bcells.ix]
# qvalsBH.df <- melt(qBH)
# colnames(qvalsBH.df) <- c("B-cell", "parameter", "qvalue")
# qvalsBH.df.ranked <- qvalsBH.df[order(qvalsBH.df$qvalue, decreasing = F),]
# qvalsBH.df.ranked[qvalsBH.df.ranked$qvalue <= 0.3,]

# Local FDR
require(fdrtool, quietly = T)
invisible(suppressMessages(fdrobj <- fdrtool(as.numeric(pvals.matrix), statistic = "pvalue", plot = F, verbose = F)))
qvals.matrix <- matrix(fdrobj$q, nrow = nrow(pvals.matrix), ncol = ncol(pvals.matrix), byrow = F)
rownames(qvals.matrix) <- rownames(pvals.matrix)
colnames(qvals.matrix) <- colnames(pvals.matrix)
rownames(qvals.matrix) <- names(df)[bcells.ix]

pvals.df <- melt(pvals.matrix)
qvals.df <- melt(qvals.matrix)
colnames(qvals.df) <- colnames(pvals.df) <- c("B-cell", "parameter", "qvalue")
qvals.df.short <- qvals.df[pvals.df$qvalue <= 0.05,]
shortlist <- qvals.df.short[order(qvals.df.short$qvalue),]



kable(shortlist, 
      digits = 3,
      row.names = T,
      caption = "Significant results") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  font_size = 10) %>%
    scroll_box(width = "100%")

Significant results

	B-cell	parameter	qvalue
87	CD19+27+IgD-38++IgG ASC	LowIntermacs	0.104
32	lymph	Sex	0.314
2	num lymph	AgeGreater60	0.348
3	lymph	AgeGreater60	0.417
45	CD27+IgD+IgM+ nonswitched memory	Sex	0.675
101	CD27+IgD+ unswitched memory	RVAD	0.721
29	CD19+27+IgD-38++IgG ASC	AgeGreater60	0.725
20	CD19CD24hiCD38-memory	AgeGreater60	0.732
61	lymph	LowIntermacs	0.732
55	CD19+CD27+CD24hi	Sex	0.738

We plotted the average across time for each of the B-cells that showed a statistically significant effect across various factors in the above one-way anova. We drew attention to specific features that induced the positive test result, by listing the model parameters with \(p<0.05\) in the multivariate fit.

require(stringr, quietly = T)
siggroups <- sapply(str_split(shortlist$parameter, ":"), function(x) x[1])
for(ii in 1:nrow(shortlist)){
    this_group <-siggroups[ii]
    this_bcell <- as.character(shortlist$`B-cell`[ii])
    cat("  \n###", as.character(shortlist$`B-cell`[ii]), "\n")
    
    sumtable <- suppressMessages(summary(models.b[[this_group]]$model[[this_bcell]]$model))
    sumtable <- as.data.frame(sumtable$coefficients)[-1, ,drop=F] # drop intercept
    if(!("Pr(>|t|)" %in% colnames(sumtable))) next()
    sigsum <- sumtable[sumtable[,"Pr(>|t|)"] <= 0.05, , drop = F]
    
    print(kable(sigsum[order(sigsum[,"Pr(>|t|)"]),,drop=F], row.names = T) %>%
              kable_styling(bootstrap_options = c("striped", 
                                                  "hover", 
                                                  "condensed",
                                                  "responsive"),
                            font_size = 12)
    )
    cat("  \n")
    
    suppressWarnings(print(plots.ts[[shortlist$`B-cell`[ii]]][[which(groups == this_group)]]))
    cat("  \n")
}

CD19+27+IgD-38++IgG ASC

	Estimate	Std. Error	df	t value	Pr(>|t|)
LowIntermacsHigh	2.078718	0.60011	17.90374	3.463894	0.0027878

lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)
SexMale	-18.44662	6.370603	17.92503	-2.895584	0.0096664

num lymph

lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)
AgeGreater60older	-15.22801	5.790664	17.53477	-2.629751	0.0172536

CD27+IgD+IgM+ nonswitched memory

	Estimate	Std. Error	df	t value	Pr(>|t|)
SexMale	-16.26667	7.449772	17.76457	-2.183512	0.0426655

CD27+IgD+ unswitched memory

	Estimate	Std. Error	df	t value	Pr(>|t|)

CD19+27+IgD-38++IgG ASC

	Estimate	Std. Error	df	t value	Pr(>|t|)

CD19CD24hiCD38-memory

	Estimate	Std. Error	df	t value	Pr(>|t|)

lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)

CD19+CD27+CD24hi

	Estimate	Std. Error	df	t value	Pr(>|t|)

Linear mixed-effect model

We analyzed the differences in B-cell levels for various features using a linear mixed effect model, with time as a continuous variable. Here we report variables that had a statistically significant variance (\(p<0.05\)) across groups, or groups at each timepoint.

suppressMessages(require(lmerTest, quietly = TRUE))
suppressMessages(require(car, quietly = TRUE))
require(reshape2, quietly = TRUE)

df.lmer <- df.HMII
names(df.lmer) <- make.names(names(df.lmer), unique = TRUE)

groupvars.ix <- c(3,4,5,7,8,9,11)
groupvars <- names(df.lmer)[groupvars.ix]

bcells.ix <- c(14:42)
bcells <- names(df.lmer)[bcells.ix]

models.b <- mclapply(groupvars, function(this_groupvar){
    models.bcells <- lapply(bcells, function(this_bcell){
        this_formula <- as.formula(paste0(this_bcell, " ~ ", this_groupvar, 
                                          " * Time + (1|PatientID)"))
        suppressMessages(suppressWarnings(this_model <- lmer(this_formula, data = droplevels(df.lmer))))
        this_anova <- Anova(this_model, type = 2)
        this_pvalues <- this_anova$`Pr(>Chisq)`
        names(this_pvalues) <- rownames(this_anova)
        #return(this_pvalues)
        return(list(model = this_model,
                    pvals = this_pvalues))
    })
    names(models.bcells) <- colnames(df)[14:42]
    pvalues <- do.call(rbind, lapply(models.bcells, function(x) x$pvals))
    rownames(pvalues) <- bcells
    #return(pvalues)
    return(list(model = models.bcells,
                pvals = pvalues))
}, mc.cores = detectCores()-1)
names(models.b) <- groupvars
pvals <- lapply(models.b, function(x) x$pvals)
# something wrong here
names(pvals) <- groupvars
pvals.matrix <- do.call(cbind, lapply(pvals, function(this_pval) this_pval[,c(1,3)]))



# Benjamini Hochberg
# qBH <- matrix(p.adjust(as.numeric(pvals.matrix), 
#                        method = "BH"), 
#               nrow = nrow(pvals.matrix), 
#               ncol = ncol(pvals.matrix), 
#               byrow = F) 
# rownames(qBH) <- rownames(pvals.matrix)
# colnames(qBH) <- colnames(pvals.matrix)
# rownames(qBH) <- names(df)[bcells.ix]
# qvalsBH.df <- melt(qBH)
# colnames(qvalsBH.df) <- c("B-cell", "parameter", "qvalue")
# qvalsBH.df.ranked <- qvalsBH.df[order(qvalsBH.df$qvalue, decreasing = F),]
# qvalsBH.df.ranked[qvalsBH.df.ranked$qvalue <= 0.3,]

# Local FDR
require(fdrtool, quietly = T)
invisible(suppressMessages(fdrobj <- fdrtool(as.numeric(pvals.matrix), statistic = "pvalue", plot = F, verbose = F)))
qvals.matrix <- matrix(fdrobj$q, nrow = nrow(pvals.matrix), ncol = ncol(pvals.matrix), byrow = F)
rownames(qvals.matrix) <- rownames(pvals.matrix)
colnames(qvals.matrix) <- colnames(pvals.matrix)
rownames(qvals.matrix) <- names(df)[bcells.ix]

pvals.df <- melt(pvals.matrix)
qvals.df <- melt(qvals.matrix)
colnames(qvals.df) <- colnames(pvals.df) <- c("B-cell", "parameter", "qvalue")
qvals.df.short <- qvals.df[pvals.df$qvalue <= 0.05,]
shortlist <- qvals.df.short[order(qvals.df.short$qvalue),]



kable(shortlist, 
      digits = 3,
      row.names = T,
      caption = "Significant results") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  font_size = 10) %>%
    scroll_box(width = "100%")

Significant results

	B-cell	parameter	qvalue
31	num lymph	AgeGreater60:Time	0.041
145	CD19+27+IgD-38++IgG ASC	LowIntermacs	0.041
2	num lymph	AgeGreater60	0.056
347	CD19+CD5+CD11b+	VAD.Indication:Time	0.060
105	CD19+24dim38dim naive mature	Sex:Time	0.217
61	lymph	Sex	0.247
344	CD19+CD5+	VAD.Indication:Time	0.319
155	CD27-38++ transitional	LowIntermacs:Time	0.331
343	CD19+CD11b+	VAD.Indication:Time	0.358
3	lymph	AgeGreater60	0.401
346	CD19+CD5+CD24hi	VAD.Indication:Time	0.413
53	CD19+CD11b+	AgeGreater60:Time	0.439
30	num Total PBMC	AgeGreater60:Time	0.508
39	CD27-38++ transitional	AgeGreater60:Time	0.537
382	CD3 of live lymph	Survival:Time	0.549
74	CD27+IgD+IgM+ nonswitched memory	Sex	0.549
92	CD3 of live lymph	Sex:Time	0.569
97	CD27-38++ transitional	Sex:Time	0.572
188	CD27+IgD+ unswitched memory	RVAD	0.574
119	lymph	LowIntermacs	0.594
57	CD19+CD5+CD11b+	AgeGreater60:Time	0.605
20	CD19CD24hiCD38-memory	AgeGreater60	0.613
29	CD19+27+IgD-38++IgG ASC	AgeGreater60	0.615
328	CD27+38++plasma blasts	VAD.Indication:Time	0.616
84	CD19+CD27+CD24hi	Sex	0.618
321	num lymph	VAD.Indication:Time	0.622
52	CD268 of +27-38++transitional	AgeGreater60:Time	0.626
399	CD19+CD268+	Survival:Time	0.627

We plotted the average across time for each of the B-cells that showed a statistically significant effect across various factors in the above mixed effect models. We drew attention to specific features that induced the positive test result, by listing the model parameters with \(p<0.05\) in the multivariate fit.

require(stringr, quietly = T)
siggroups <- sapply(str_split(shortlist$parameter, ":"), function(x) x[1])
for(ii in 1:nrow(shortlist)){
    this_group <-siggroups[ii]
    this_bcell <- as.character(shortlist$`B-cell`[ii])
    cat("  \n###", as.character(shortlist$`B-cell`[ii]), "\n")
    
    sumtable <- suppressMessages(summary(models.b[[this_group]]$model[[this_bcell]]$model))
    sumtable <- as.data.frame(sumtable$coefficients)[-1, ,drop=F] # drop intercept
    if(!("Pr(>|t|)" %in% colnames(sumtable))) next()
    sigsum <- sumtable[sumtable[,"Pr(>|t|)"] <= 0.05, , drop = F]
    
    print(kable(sigsum[order(sigsum[,"Pr(>|t|)"]),,drop=F], row.names = T) %>%
              kable_styling(bootstrap_options = c("striped", 
                                                  "hover", 
                                                  "condensed",
                                                  "responsive"),
                            font_size = 12)
    )
    cat("  \n")
    
    suppressWarnings(print(plots.ts[[shortlist$`B-cell`[ii]]][[which(groups == this_group)]]))
    cat("  \n")
}

num lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	8206.115	1791.316	1287.688	4.581054	0.0000051
AgeGreater60older:Time	-8800.327	2381.169	2274.895	-3.695802	0.0002244

CD19+27+IgD-38++IgG ASC

	Estimate	Std. Error	df	t value	Pr(>|t|)
LowIntermacsHigh	3.058696	0.774233	45.35361	3.950614	0.0002697

num lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	8206.115	1791.316	1287.688	4.581054	0.0000051
AgeGreater60older:Time	-8800.327	2381.169	2274.895	-3.695802	0.0002244

CD19+CD5+CD11b+

	Estimate	Std. Error	df	t value	Pr(>|t|)
VAD.IndicationDT:Time	0.7535101	0.2208725	66.53183	3.411516	0.0011029

CD19+24dim38dim naive mature

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	-0.8161430	0.2628953	64.83737	-3.104441	0.0028253
SexMale:Time	0.9045098	0.3076620	64.86978	2.939946	0.0045445

lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)
SexMale	-17.61531	7.156886	25.50523	-2.46131	0.0209393

CD19+CD5+

	Estimate	Std. Error	df	t value	Pr(>|t|)
VAD.IndicationDT:Time	0.7495301	0.2761314	67.17233	2.714397	0.0084299

CD27-38++ transitional

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	0.2728189	0.0943080	68.33055	2.892849	0.0051174
LowIntermacsHigh:Time	-0.3029055	0.1126997	67.21515	-2.687721	0.0090598

CD19+CD11b+

	Estimate	Std. Error	df	t value	Pr(>|t|)
VAD.IndicationDT:Time	0.5911116	0.2251629	65.39624	2.625262	0.0107693

lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	0.7428399	0.3610038	69.6636	2.057706	0.0433624

CD19+CD5+CD24hi

	Estimate	Std. Error	df	t value	Pr(>|t|)
VAD.IndicationDT:Time	0.4514841	0.1801168	68.30036	2.506618	0.0145744

CD19+CD11b+

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	0.5499941	0.1571744	65.94444	3.499260	0.0008412
AgeGreater60older:Time	-0.5051776	0.2059164	65.50982	-2.453314	0.0168239

num Total PBMC

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	8539.807	3421.824	8916.676	2.495688	0.0125892
AgeGreater60older:Time	-10420.156	4502.010	15842.034	-2.314556	0.0206499

CD27-38++ transitional

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	0.1924463	0.0792502	69.27151	2.428340	0.0177687
AgeGreater60older:Time	-0.2343863	0.1043092	67.73025	-2.247034	0.0279003

CD3 of live lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)
Survivaldead:Time	-1.829251	0.8242817	66.32952	-2.219206	0.0298979

CD27+IgD+IgM+ nonswitched memory

	Estimate	Std. Error	df	t value	Pr(>|t|)
SexMale	-20.7260608	7.6448089	21.84155	-2.711129	0.0128054
Time	-0.9006538	0.3895893	66.29131	-2.311803	0.0239045

CD3 of live lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)
SexMale:Time	-1.143545	0.5277063	66.36504	-2.167011	0.0338279

CD27-38++ transitional

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	0.2418057	0.0994538	65.99621	2.431336	0.0177663
SexMale:Time	-0.2512294	0.1163638	66.05080	-2.159000	0.0344873

CD27+IgD+ unswitched memory

	Estimate	Std. Error	df	t value	Pr(>|t|)
RVADYes	8.138537	3.426992	21.41969	2.374834	0.0269627

lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)

CD19+CD5+CD11b+

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	0.5382593	0.1602283	67.49846	3.359327	0.0012883
AgeGreater60older:Time	-0.4334330	0.2101801	66.78306	-2.062198	0.0430827

CD19CD24hiCD38-memory

	Estimate	Std. Error	df	t value	Pr(>|t|)

CD19+27+IgD-38++IgG ASC

	Estimate	Std. Error	df	t value	Pr(>|t|)
AgeGreater60older	2.044292	0.8970807	39.13157	2.278828	0.0282137

CD27+38++plasma blasts

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	0.1018402	0.0463711	65.98620	2.196200	0.0315970
VAD.IndicationDT:Time	-0.1781890	0.0880343	66.32797	-2.024086	0.0469921

CD19+CD27+CD24hi

	Estimate	Std. Error	df	t value	Pr(>|t|)

num lymph

CD268 of +27-38++transitional

	Estimate	Std. Error	df	t value	Pr(>|t|)

CD19+CD268+

	Estimate	Std. Error	df	t value	Pr(>|t|)
Time	-0.5331745	0.2039893	64.53338	-2.613738	0.0111329

Two-way repeated measures anova

We analyzed the differences in B-cell levels for various features using a two-way repeated measures ANOVA. Here we report variables that had a statistically significant variance (\(p<0.05\)) across groups, or groups at each timepoint.

suppressMessages(require(lmerTest, quietly = TRUE))
suppressMessages(require(car, quietly = TRUE))
require(reshape2, quietly = TRUE)

df.lmer <- df.HMII
names(df.lmer) <- make.names(names(df.lmer), unique = TRUE)

groupvars.ix <- c(3,4,5,7,8,9,11)
groupvars <- names(df.lmer)[groupvars.ix]

bcells.ix <- c(14:42)
bcells <- names(df.lmer)[bcells.ix]

models.b <- mclapply(groupvars, function(this_groupvar){
    models.bcells <- lapply(bcells, function(this_bcell){
        this_formula <- as.formula(paste0(this_bcell, " ~ ", this_groupvar, 
                                          " * factor(Time) + (1|PatientID)"))
        suppressMessages(suppressWarnings(this_model <- lmer(this_formula, data = droplevels(df.lmer))))
        this_anova <- Anova(this_model, type = 2)
        this_pvalues <- this_anova$`Pr(>Chisq)`
        names(this_pvalues) <- rownames(this_anova)
        #return(this_pvalues)
        return(list(model = this_model,
                    pvals = this_pvalues))
    })
    names(models.bcells) <- colnames(df)[14:42]
    pvalues <- do.call(rbind, lapply(models.bcells, function(x) x$pvals))
    rownames(pvalues) <- bcells
    #return(pvalues)
    return(list(model = models.bcells,
                pvals = pvalues))
}, mc.cores = detectCores()-1)
names(models.b) <- groupvars
pvals <- lapply(models.b, function(x) x$pvals)
# something wrong here
names(pvals) <- groupvars
pvals.matrix <- do.call(cbind, lapply(pvals, function(this_pval) this_pval[,c(1,3)]))



# Benjamini Hochberg
# qBH <- matrix(p.adjust(as.numeric(pvals.matrix), 
#                        method = "BH"), 
#               nrow = nrow(pvals.matrix), 
#               ncol = ncol(pvals.matrix), 
#               byrow = F) 
# rownames(qBH) <- rownames(pvals.matrix)
# colnames(qBH) <- colnames(pvals.matrix)
# rownames(qBH) <- names(df)[bcells.ix]
# qvalsBH.df <- melt(qBH)
# colnames(qvalsBH.df) <- c("B-cell", "parameter", "qvalue")
# qvalsBH.df.ranked <- qvalsBH.df[order(qvalsBH.df$qvalue, decreasing = F),]
# qvalsBH.df.ranked[qvalsBH.df.ranked$qvalue <= 0.3,]

# Local FDR
require(fdrtool, quietly = T)
invisible(suppressMessages(fdrobj <- fdrtool(as.numeric(pvals.matrix), statistic = "pvalue", plot = F, verbose = F)))
qvals.matrix <- matrix(fdrobj$q, nrow = nrow(pvals.matrix), ncol = ncol(pvals.matrix), byrow = F)
rownames(qvals.matrix) <- rownames(pvals.matrix)
colnames(qvals.matrix) <- colnames(pvals.matrix)
rownames(qvals.matrix) <- names(df)[bcells.ix]

pvals.df <- melt(pvals.matrix)
qvals.df <- melt(qvals.matrix)
colnames(qvals.df) <- colnames(pvals.df) <- c("B-cell", "parameter", "qvalue")
qvals.df.short <- qvals.df[pvals.df$qvalue <= 0.05,]
shortlist <- qvals.df.short[order(qvals.df.short$qvalue),]



kable(shortlist, 
      digits = 3,
      row.names = T,
      caption = "Significant results") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  font_size = 10) %>%
    scroll_box(width = "100%")

Significant results

	B-cell	parameter	qvalue
2	num lymph	AgeGreater60	0.040
145	CD19+27+IgD-38++IgG ASC	LowIntermacs	0.040
52	CD268 of +27-38++transitional	AgeGreater60:factor(Time)	0.073
276	CD27+IgD-IgM+ switched memory	Sensitized:factor(Time)	0.081
61	lymph	Sex	0.200
348	CD19+27+IgD-38++IgG ASC	VAD.Indication:factor(Time)	0.213
232	CD19+27+IgD-38++IgG ASC	RVAD:factor(Time)	0.230
155	CD27-38++ transitional	LowIntermacs:factor(Time)	0.254
3	lymph	AgeGreater60	0.257
347	CD19+CD5+CD11b+	VAD.Indication:factor(Time)	0.377
31	num lymph	AgeGreater60:factor(Time)	0.378
101	CD27+IgD+ unswitched memory	Sex:factor(Time)	0.417
40	CD27-IgD+ mature naive	AgeGreater60:factor(Time)	0.438
119	lymph	LowIntermacs	0.528
321	num lymph	VAD.Indication:factor(Time)	0.592
74	CD27+IgD+IgM+ nonswitched memory	Sex	0.616
217	CD27+IgD+ unswitched memory	RVAD:factor(Time)	0.621
103	CD27+IgD+IgM+ nonswitched memory	Sex:factor(Time)	0.632
225	CD19+CD268+	RVAD:factor(Time)	0.643
399	CD19+CD268+	Survival:factor(Time)	0.675
344	CD19+CD5+	VAD.Indication:factor(Time)	0.677
224	CD19+27-38+CD5+transitionals	RVAD:factor(Time)	0.705
51	CD19+CD268+	AgeGreater60:factor(Time)	0.705
188	CD27+IgD+ unswitched memory	RVAD	0.708
29	CD19+27+IgD-38++IgG ASC	AgeGreater60	0.727
20	CD19CD24hiCD38-memory	AgeGreater60	0.732

We plotted the average across time for each of the B-cells that showed a statistically significant effect across various factors in the above mixed effect models. We drew attention to specific features that induced the positive test result, by listing the model parameters with \(p<0.05\) in the multivariate fit.

require(stringr, quietly = T)
siggroups <- sapply(str_split(shortlist$parameter, ":"), function(x) x[1])
for(ii in 1:nrow(shortlist)){
    this_group <-siggroups[ii]
    this_bcell <- as.character(shortlist$`B-cell`[ii])
    cat("  \n###", as.character(shortlist$`B-cell`[ii]), "\n")
    
    sumtable <- suppressMessages(summary(models.b[[this_group]]$model[[this_bcell]]$model))
    sumtable <- as.data.frame(sumtable$coefficients)[-1, ,drop=F] # drop intercept
    if(!("Pr(>|t|)" %in% colnames(sumtable))) next()
    sigsum <- sumtable[sumtable[,"Pr(>|t|)"] <= 0.05, , drop = F]
    
    print(kable(sigsum[order(sigsum[,"Pr(>|t|)"]),,drop=F], row.names = T) %>%
              kable_styling(bootstrap_options = c("striped", 
                                                  "hover", 
                                                  "condensed",
                                                  "responsive"),
                            font_size = 12)
    )
    cat("  \n")
    
    suppressWarnings(print(plots.ts[[shortlist$`B-cell`[ii]]][[which(groups == this_group)]]))
    cat("  \n")
}

num lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)
factor(Time)21	169798.9	44963.12	1780.755	3.776403	0.0001643
AgeGreater60older:factor(Time)21	-214015.6	62340.16	1651.364	-3.433029	0.0006116
AgeGreater60older:factor(Time)5	-136885.7	61168.57	2679.738	-2.237844	0.0253131
factor(Time)14	128214.1	57636.41	1470.433	2.224532	0.0262639

CD19+27+IgD-38++IgG ASC

	Estimate	Std. Error	df	t value	Pr(>|t|)
LowIntermacsHigh	3.095	1.279721	69	2.418497	0.0182277

CD268 of +27-38++transitional

	Estimate	Std. Error	df	t value	Pr(>|t|)
AgeGreater60older:factor(Time)8	-27.22177	10.33554	54.58048	-2.633803	0.0109633
factor(Time)14	-24.92941	12.43679	55.44507	-2.004489	0.0499111

CD27+IgD-IgM+ switched memory

	Estimate	Std. Error	df	t value	Pr(>|t|)
SensitizedYes:factor(Time)3	17.965877	5.302290	23.31974	3.388324	0.0024949
factor(Time)5	9.747852	4.547781	23.30369	2.143431	0.0427391

lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)
SexMale	-22.22807	9.001586	53.16972	-2.46935	0.0167811

CD19+27+IgD-38++IgG ASC

	Estimate	Std. Error	df	t value	Pr(>|t|)
factor(Time)5	4.722429	1.330913	60.14685	3.548263	0.0007592
VAD.IndicationDT:factor(Time)5	-6.577445	1.874482	57.09611	-3.508940	0.0008850

CD19+27+IgD-38++IgG ASC

	Estimate	Std. Error	df	t value	Pr(>|t|)
RVADYes:factor(Time)5	11.55470	2.912778	59.01744	3.966900	0.0001997
factor(Time)14	-2.17979	1.043488	56.60858	-2.088945	0.0412202

CD27-38++ transitional

	Estimate	Std. Error	df	t value	Pr(>|t|)
LowIntermacsHigh:factor(Time)21	-11.380113	2.894440	59.65377	-3.931715	0.0002222
factor(Time)21	9.440964	2.450218	61.33982	3.853112	0.0002815

lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)

CD19+CD5+CD11b+

	Estimate	Std. Error	df	t value	Pr(>|t|)
VAD.IndicationDT:factor(Time)21	18.60843	5.699112	56.68547	3.265145	0.0018584
VAD.IndicationDT:factor(Time)8	10.32734	4.156816	55.37116	2.484436	0.0160275

num lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)
factor(Time)21	169798.9	44963.12	1780.755	3.776403	0.0001643
AgeGreater60older:factor(Time)21	-214015.6	62340.16	1651.364	-3.433029	0.0006116
AgeGreater60older:factor(Time)5	-136885.7	61168.57	2679.738	-2.237844	0.0253131
factor(Time)14	128214.1	57636.41	1470.433	2.224532	0.0262639

CD27+IgD+ unswitched memory

	Estimate	Std. Error	df	t value	Pr(>|t|)
SexMale:factor(Time)8	9.848232	2.746704	54.41204	3.585472	0.0007200
factor(Time)8	-8.159681	2.352215	54.37377	-3.468936	0.0010300
factor(Time)14	-8.101354	3.062814	55.12449	-2.645069	0.0106199
SexMale	-10.353852	3.922084	26.14988	-2.639885	0.0138005
factor(Time)5	-7.781354	3.062814	55.12449	-2.540590	0.0139123
SexMale:factor(Time)5	8.682914	3.494502	55.06965	2.484736	0.0160327
SexMale:factor(Time)3	5.898797	2.849382	54.68689	2.070202	0.0431688
SexMale:factor(Time)14	7.448953	3.625289	54.98976	2.054720	0.0446698

CD27-IgD+ mature naive

	Estimate	Std. Error	df	t value	Pr(>|t|)
factor(Time)14	-19.34727	6.403688	54.77954	-3.021270	0.0038210
AgeGreater60older:factor(Time)14	17.50253	7.722534	54.66849	2.266424	0.0274049

lymph

	Estimate	Std. Error	df	t value	Pr(>|t|)

num lymph

CD27+IgD+IgM+ nonswitched memory

	Estimate	Std. Error	df	t value	Pr(>|t|)
SexMale	-31.50752	9.007301	37.48741	-3.497997	0.0012251
SexMale:factor(Time)8	27.02760	8.735159	54.73057	3.094116	0.0031071
SexMale:factor(Time)3	24.44017	9.046112	55.34852	2.701732	0.0091386
factor(Time)8	-20.00903	7.482252	54.65698	-2.674199	0.0098595
factor(Time)3	-15.59653	7.482252	54.65698	-2.084470	0.0418040

CD27+IgD+ unswitched memory

	Estimate	Std. Error	df	t value	Pr(>|t|)
RVADYes:factor(Time)8	-10.16212	3.121996	55.38516	-3.255008	0.0019363
RVADYes	11.57901	3.999725	34.08235	2.894952	0.0065730

CD27+IgD+IgM+ nonswitched memory

	Estimate	Std. Error	df	t value	Pr(>|t|)
SexMale	-31.50752	9.007301	37.48741	-3.497997	0.0012251
SexMale:factor(Time)8	27.02760	8.735159	54.73057	3.094116	0.0031071
SexMale:factor(Time)3	24.44017	9.046112	55.34852	2.701732	0.0091386
factor(Time)8	-20.00903	7.482252	54.65698	-2.674199	0.0098595
factor(Time)3	-15.59653	7.482252	54.65698	-2.084470	0.0418040

CD19+CD268+

	Estimate	Std. Error	df	t value	Pr(>|t|)
RVADYes:factor(Time)3	-27.99543	10.480918	56.66602	-2.671086	0.0098514
RVADYes:factor(Time)8	-20.33650	9.017817	55.39973	-2.255146	0.0280950

CD19+CD268+

	Estimate	Std. Error	df	t value	Pr(>|t|)
Survivaldead:factor(Time)8	-21.73092	8.069861	55.34955	-2.692849	0.0093560
Survivaldead:factor(Time)3	-18.76413	8.701780	55.67137	-2.156355	0.0353912

CD19+CD5+

	Estimate	Std. Error	df	t value	Pr(>|t|)
VAD.IndicationDT:factor(Time)8	13.77699	5.156496	55.63351	2.671773	0.0098785
VAD.IndicationDT:factor(Time)21	18.57915	7.056032	57.34507	2.633087	0.0108559

CD19+27-38+CD5+transitionals

	Estimate	Std. Error	df	t value	Pr(>|t|)

CD19+CD268+

	Estimate	Std. Error	df	t value	Pr(>|t|)
factor(Time)14	-34.03871	8.813153	54.50989	-3.862263	0.0003004
AgeGreater60older:factor(Time)14	34.33335	10.627446	54.41429	3.230631	0.0020977
factor(Time)21	-17.19694	6.912037	54.73282	-2.487970	0.0159221
AgeGreater60older:factor(Time)3	19.32183	8.006711	54.60648	2.413204	0.0191971
factor(Time)8	-11.89000	5.511085	53.89561	-2.157470	0.0354450

CD27+IgD+ unswitched memory

	Estimate	Std. Error	df	t value	Pr(>|t|)
RVADYes:factor(Time)8	-10.16212	3.121996	55.38516	-3.255008	0.0019363
RVADYes	11.57901	3.999725	34.08235	2.894952	0.0065730

CD19+27+IgD-38++IgG ASC

	Estimate	Std. Error	df	t value	Pr(>|t|)

CD19CD24hiCD38-memory

	Estimate	Std. Error	df	t value	Pr(>|t|)